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THREE  PARTICLE  EFFECTS  IN  THE  PAIR 

h 
DISTRIBUTION  FUNCTION  FOR  HE  GAS 


Harry  Frederick  Jordan,  Ph.D. 

Department  of  Physics 
University  of  Illinois,  19^7 


The  diagonal  elements  of  the  3 -particle  density  matrix 

h 
for  He  have  been  formulated  in  terms  of  Wiener  functional  integrals 

in  order  to  compute  both  that  part  of  the  pair  distribution  function 

linear  in  the  density  and  the  third  virial  coefficient.   The  3-particle 

interaction  energy  was  taken  to  be  a  sum  of  Lennard -Jones  (12-6) 

potentials  between  pairs  of  particles.   Evaluation  of  the  Wiener 

integrals  was  carried  out  by  a  Monte  Carlo  sampling  technique  to  yield 

the  pair  distribution  function  and  third  virial  coefficient  for 

temperatures  from  273  K  to  5  K.   Comparison  of  the  linear  portion  of 

the  density  expansion  of  the  pair  distribution  function  gave 

qualitative  agreement  with  experimental  values  for  the  pair  distribution 

k 
function  in  liquid  He  measured  by  neutron  diffraction.   The  values 

of  the  third  virial  coefficient  calculated  by  this  method  are 

consistently  smaller  than  the  experimental  values  and  consistently 

larger  than  the  classical  values  of  this  quantity  over  the  given 

temperature  range . 


1 .   INTRODUCTION 

In  the  past  helium  gas  has  "been  a  common  substance  for  use  in 
the  comparison  with  experiment  of  the  results  of  theoretical  models  and 
techniques  in  quantum  statistical  mechanics.   At  the  time  quantum  theory 
began  to  be  well  understood  the  statistical  mechanical  theory  of  gases 
had  developed  further  than  the  theory  of  either  solids  or  liquids.   The 
use  of  helium  as  a  test  substance  arose  due  to  the  fact  that  of  the 
monatomic  gases,  for  which  the  application  of  statistical  mechanics 
presents  the  least  complexity,  helium  is  the  only  one  having  small 
enough  atomic  weight  to  exhibit  sizeable  quantum  effects  in  its 
thermodynamic  properties  at  moderate  experimental  temperatures  and 
pressures. 

The  second  virial  coefficient  B(T)  has  received  a  great  deal 
of  attention,  and  calculations  of  B(T)  have  been  made  using  several 
different  forms  of  the  intermolecular  potential  (for  a  survey  of  much 
of  this  work  see  Hirshfelder,  Curtis,  and  Bird  [l]).   The  two  methods 
most  frequently  applied  to  the  calculation  of  B(T)  are  the  expression 
of  B(T)  in  terms  of  the  phase  shifts  of  2-particle  collision  theory 
and  the  expansion  of  B(T)  about  the  classical  limit  as  a  power  series 
in  Planck's  constant  h  (the  Wigner-Kirkwood  expansion).  The  latter 
method  is  useful  only  for  high  temperatures  since  the  series  converges 
slowly,  if  at  all,  for  temperatures  less  than  about  ^0  K  in  He  . 


Few  quantum  mechanical  calculations  have  been  made  for  higher 
order  virial  coefficients.  Larsen  [2]  and  Pais  and  Uhlenbeck  [3]  have 
done  work  on  the  third  virial  coefficient.  Both  calculations  have  made 
use  of  the  "binary  collision  expansion  of  Lee  and  Yang  [h] ,    for  which 
the  convergence  properties  are  not  well  understood.   The  Wigner- 
Kirkwood  expansion  is  applicable  to  systems  of  more  than  two  particles 
and  can  thus  be  used  to  compute  higher  order  virial  coefficients,  but 
the  calculations  are  more  complex,  and  the  convergence  difficulties  are 
the  same  as  in  the  2 -particle  case.   The  phase  shift  approach,  of 
course,  is  only  valid  for  quantities  which  may  be  expressed  in  terms  of 
2-body  collision  theory. 

Recently  two  calculations  of  the  trace  of  the  2 -particle 
density  matrix,  a  2-body  property  more  general  than  the  second  virial, 
have  been  made  in  ways  which  admit  of  generalization  to  the  N-particle 
problem  for  N  >  2.   One  of  these  [5]  involves  the  numerical  calculation 
of  the  eigenf unctions  and  eigenvalues  of  the  Schro'dinger  equation  for 
the  relative  coordinate  and  explicit  evaluation  of  the  sum  over  states 
involved  in  the  density  matrix.   The  other  [6]  makes  use  of  a  path 
integral  formulation  for  the  quantum  mechanical  density  matrix  and  a 
Monte  Carlo  evaluation  of  the  multiple  integrals  which  arise  in  the 
evaluation  of  these  path  integrals.  This  thesis  concerns  the 
application  of  the  latter  method  to  the  3 -particle  problem  in  order 
to  evaluate  both  the  3_particle  term  in  the  pair  distribution  function 
and  the  third  virial  coefficient, 


2.   THEORY  OF  THE  PAIR  DISTRIBUTION  FUNCTION 

2.1   Definition  of  the  Pair  Distribution  Function  in  Terms  of  the 
Density  Matrix 

k 

To  clarify  the  mathematical  model  to  be  used  for  He  we  assume 

that  the  system  of  interest  is  a  gas  of  N  identical  particles  in  a 

volume  V  in  the  limit  as  N,  V  -+  °°  with  n  =  n/v  held  constant.   The 

particles,  the  coordinates  of  which  are  denoted  by  i,  i=l,2,...,N,  are 

assumed  to  interact  by  means  of  an  additive,  spherically  symmetric  pair 

potential  V(r.  . ) ,   where  r.  .  =  | i-jj ,  so  that  the  interaction  energy  for 
•^  J         -L  J 

the  N-particle  system  is 


U(1,2,...,N)  =  |  Z  V(r   ).  (2.1.1) 


In  the  numerical  calculations  the  potential  function  will  be  taken  to 
be  the  Lennard- Jones  (12-6)  potential 


V(r)  =  k-K. 


l  fe  r  -  (? 


(2,1.2) 


where  k  =  ik.Ok  x  10  '      erg  and  a  =   2.56  X  10    cm.  This  potential  is 
taken  because  it  is  computationally  simple  and  gives  almost  as  good  a 

fit  to  second  virial  coefficient  data,  at  least,  as  any  of  the  other 

k 

potentials  proposed  for  He  .   The  constants  k  and  a   are  taken  after 

de  Boer  and  Michels  [7]  to  give  a  good  fit  to  high  temperature  second 
virial  coefficients.  The  assumption  of  an  additive  pair  potential  is 
not  essential  to  the  theoretical  development  below.   In  order  to 


produce  results  for  more  realistic  potentials  one  need  only  replace  the 
3-partiele  potential  energy,  U  (1,2,3),  in  the  final  formulae  by 
the  appropriate  3-particle  interaction. 

The  pair  distribution  function  n  (l,2)d  1  d  2  is  defined  as 

the  probability  that  one  of  the  particles  of  the  system  is  found  in  the 

3 
volume  element  d  1  at  1  while  some  other  particle  is  simultaneously  in 

3 

volume  d  2  at  2.   In  the  gaseous  system  under  consideration,  in  fact 

in  any  homogeneous  fluid,  n  (1,2)  must  depend  only  on  the  separation  r,p 
between  the  particles  except  in  the  region  near  the  boundaries  of  the 
system  which  becomes  negligible  as  V  -»  <*.   In  a  classical  gas  of  non- 
interacting  particles  the  probability  of  finding  a  particle  in  a  given 
volume  element  becomes  a  constant  independent  of  the  positions  of  the 
other  particles,  so  in  this  case 

/.  Qx    N(N-l)  ,0  ,  -v 

no'l>V    =    ~p >  (2.1-3) 

V 


lim  n2U,2)  -  n2  .  (2.1.4) 

N,V-oa 

n  =  const , 


Quantum  mechanically  the  situation  is  complicated  by  the  exchange 
correlation  between  identical  particles,  but  since  this  effect  falls 
off  as  (see  Section  3*2)  exp[ -const,  x  (r  A)  ],  Eqs .  (2.1.3)  and 
(2.1.4)  still  hold  in  the  limit  r, p  »  X,  where 


is  the  thermal  wavelength  of  a  particle  of  mass  m  at  inverse 
temperature  £  =  l/kT. 

If  PN(1,2, . ..,N)d3l  d32  ...  d3N  is  the  probability  that 

3 

particle  i  is  in  volume  d  i  at  i  for  i  from  1  to  N,  then  the  pair 

distribution  function  is  given  by 

n2(l,2)  =  l(H)J...JPI(l)2,...)N)d33  d\    ...  d3N.    (2.1.5) 


The  conf igurational  probability  P  can  be  simply  expressed  in  terms  of 
the  diagonal  elements  of  the  density  matrix  for  the  quantum  mechanical 
canonical  ensemble 

P  (1,2,..., N)  =        <l,2,...,N|e    M|l,2,...,N>.     (2.1.6) 
^  \ 


Here  H^  is  the  Hamiltonian  operator  of  the  N-particle  system 

N  £ 


in  Jd 
HN=  ~Z  f^  +  Vi^---*3)  >  (2-1^) 


Z   is  the  partition  function 


Z  =  -~=  /  ...  /  <a,2,...,N|e   ^|1,2,...,N>  d3l  d32  .  .  .  d3N  , 
NX"3  J      V  J 


(2.1.8) 


and  the  density  matrix  is  given  in  terms  of  the  energy  eigenf unctions 
d)  (1,2,  .  .  .  ,N)  of  the  system  as 


<1' ,2* ,...,$*  \e  |l,2,...,N> 

_Qtp 

=  NI\3NZ  M(V,  2  ',...,  N')(j)  ,(1,2,...,  N)e   J   ,  (2.1.9) 

J  J  J 

where  H„(t> .  =  E.d)..   The  pair  distribution  function  can  then  be  written 


n2(i'S)  =  " riN~  J  '*'/<aJ'2,...,H|e   N|1,2,...,N> 

d        (N-2)!^ZN  J   v  J 


x  d3^  d3i+  ...d3N  .  (2.1.10) 


2.2   The  Effect  of  Spin  and  Statistics  on  the  Density  Matrix 

It  should  be  noted  that  for  a  system  of  identical  particles 
the  sum  in  Eq.  (2.1.9)  extends  only  over  properly  symmetrized  eigen- 
states.   It  will  be  useful  in  the  following  to  express  density  matrix 
elements  for  Bose  and  Fermi  statistics  in  terms  of  matrix  elements  for 
Boltzmann  statistics  and  thus  to  deal  with  unrestricted  sums  over 
eigenstates.   If  we  let  the  subscripts  +  denote  Bose  statistics  and 
-  denote  Fermi  statistics,  then  the  symmetrized  N-particle  states  for 
both  cases  can  be  written 


+ 


,2,...,N>  =  Zp  — r  |P(1,2,...,N)>  ,         (2.2.1) 
-    *  Jul 


where  P  is  one  of  the  NJ  permutations  of  N  elements  and  ap  is  defined 

by 


+     , 


(J-r.   = 


+1  ,  if  P  is  an  even  permutation,  (2.2.2) 


P  "    \ 

-1  ,  if  P  is  an  odd  permutation. 

The  Bose  and  Fermi  matrix  elements  can  thus  "be  written  in  terms  of  the 
Boltzmann  matrix  elements  as 


<l«,2',...,N'|e   iN|l,2,...,N>+ 


=  ZpI£pIap„apl<P"(l,,2',...,N,)|e    |P' (1,2, . . . ,N)>, 


Zpap  <P(l',2',...,N')|e     | 1,2, . . . ,N>  ,  (2.2.3) 


since  Z  ;iZ  ,  =  NiZ  ,  where  P  =  P"P'  .  When  space  averages  are  taken 
over  the  Boltzmann  matrix  elements  on  the  right  of  Eq.  (2.2.3)  only 
1/N!  of  phase  space  should  be  taken.   To  allow  the  integration  to 
extend  over  all  space,  we  define  the  Boltzmann  matrix  elements  without 
the  factor  of  NJ  appearing  in  the  corresponding  Bose  and  Fermi 
quantities, 


<1',2',...,N'  |e     |l,2,...,N> 


=  \^E   ty*AV,2\...,W)tyAl,2,...,K)e        J.  (2.2.4) 

all  j  J  J 

k 

Although  the  above  expressions  are  correct  for  He  ,  which  has 

spin  zero,  modifications  are  needed  for  particles  having  nonzero  spin. 

If  the  particles  each  have  spin  S,  then  an  N-particle  state, 

| l,s  ,2,s  , . . . ,N,s  >,  should  specify  a  spin  coordinate,  s.,  as  well  as 

a  space  coordinate,  i,  for  each  particle.   Since  the  Hamiltonian  is 

spin  independent,  matrix  elements  between  spin  states  yield  Kronecker 

deltas, 


<L,,s;[,2,,s£,...,]J,,s£|e    |l,S;L,2,s2,...,N,sN> 


-f3HN 
-  S(s{,s1)...6(s]^,sN)<l',£V..,N'|e     |l,2,  .  . .  ,N>  .       (2.2.5) 


If  one  then  averages  over  spin  in  Eq.  (2.2.3),  the  spin  averaged 
diagonal  density  matrix  element  becomes 


<i>2,...,N|e   S|1,2,...,H>  = ^EZ  ...  ZE  a  6(s  ,Ps  )...5(VPs^) 

-   (2S+1)   12    SN 


X  <P(1,2, ...,N)|e     |l,2, ...,N>  .  (2.2.6) 


Noting  that 


Z  5(s.,Ps.)  =  2S+1,   if  s.  =  Ps.  , 
s.    i   1  i     i 


=  1,  otherwise  , 


one  finds  that  if  the  permutation  P  is  decomposed  into  cycles,  the 

"PHN, 
term  in  <P(l,2, . . . ,N) |e     |l,2, ...,N>  contains  a  factor  of 

(2S+1)       for  each  cycle  of  length  i.   As  a  specific  example  the 

diagonal  element  of  the  3-particle  density  matrix  is 


<l,2,3|e   ^|l,2,p  =  <l,2,3|e   ^|l,2,p 


1  _PH^  1  _PH* 


1  ~^H^  1  ~^H^ 

±  ^TI  <^S,l|e    3|1,2,5>  +  2  <2,^,l|e    J|l,2,p 

V2S+1) 


1  "PH^ 

+  P  <3»i>2|e    3|1,2,3>  .  (2.2.7) 

(2S+1) 


2 .3  Cluster  Expansion  for  the  Pair  Distribution  Function 

De  Boer  [8]  has  shown  how  the  quantum  mechanical  cluster 
expansion  of  Kahn  and  Uhlenbeck  [9]  can  he  used  to  express  the  pair 
distribution  function  for  a  gas  in  terms  of  a  power  series  in  the 
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density.   It  turns  out  that  the  term  in  the  j-th  power  of  the  density 
depends  only  on  density  matrix  elements  for  systems  of  j+2  or  fewer 
particles.   The  following  brief  summary  of  definitions  and  results  draws 
from  the  treatment  in  de  Boer's  review  article  [8]. 

Denote  the  diagonal  elements  of  the  density  matrix  by 


W  (1,2,...,N)  =  <l,2,...,N|e   N|l,2,...,N>  .      (2.3-1) 


Let  r    stand  for  the  set  (1,2, ...,h)  of  coordinates  of  a  group  of  h 
particles  and  define  Ursell  functions  U«  and  modified  Ursell  functions 
\j\       by  the  relations 


u1(i)   -  w1(i)   , 

u2(i,4)  =  w2(i,j)  -  w1(i)v1(Jl)  , 

u3<  i,j,k)  =  w3(i,Aj,k)  -  w2(i,Ji)w1(k)  -  w2(i,k)w1(Ji) 

-  w2(4,^)W1(i)  +  2W]_  i)W1(Ji)W1(k)  ,  etc.  (2.3-2) 
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u(h)(r(h))=W<r(h)) 

1   v~       tr~     , 

U<h)(r(h),l)  =  Wh+1(r(h),l)  -Wh(r(h')Wl(l)  _ 

-  "2(i.j)Vs(h))  +  2V£(h))wi(i)V<i)  <  etc-  (2-3-3) 


The  Ursell  functions  have  the  property  that  if  the  particles  of  the 
system  are  separated  into  two  groups  such  that  the  minimum  distance 
between  the  groups  is  much  greater  than  both  the  thermal  wavelength  and 
the  range  of  the  interparticle  potential,  then  any  Ursell  function 
containing  as  arguments  coordinates  of  particles  in  each  group  vanishes, 
The  modified  Ursell  functions  have  the  same  property  provided  that  the 
set  of  h  particles  specified  by  r    is  considered  as  a  single  particle 
in  the  separation  into  groups.   The  Ursell  functions  are  useful  in 
treating  a  gas  since  in  this  case  the  average  interparticle  spacing  is 
large  enough  that  only  those  functions  with  a  small  number  of  arguments 
differ  from  zero. 

The  definition  of  the  modified  Ursell  functions  leads  to  an 
expression  for  W  in  terms  of  these  functions 
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Vi>2,  — ,N)   =  U[h)(£(h))WN-h(^l^'^N) 


v        (h)/  (h)      .  T  , 

f    u2  fef  ^i)\_h-i^h+i--'£i-i^i+i---%) 


N 

i=h+l 


+       E       UI.    \r/    ,r.,r.)W.T  ,'(r,    .,..., r.    .  ,r.    ,,..., r.    ,,r.   ......  r„) 

•j-_h   -i    3  ~i  ~j     N-h-2v~h+l'        7~i-l'~i+l/        '~j-1'~j+1'        ,~N/ 


+    ...  (2.3.4) 


where  r.  is  used  for  the  coordinate  of  the  i-th  particle.   If  we  define 
~i 

cluster  integrals  tA   by 

bto(PoO)  _  i  r    ru(h)(r(n)  r       r     j 

Di     -~      '      £IJ        \r        J  H  l~      '~h+l,",,~h+i-i; 


X  d3r,  nd3r,  0  ...  d3ru  .  .  ,  (2.3*5) 

~h+l  ~h+2       ~h+f-l  *  N  "*  " 


then  using  Eq.  (2.3»*0  with  h  -   2  the  pair  distribution  function  can 
be  expressed  as 
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n2(l,2)  =  ^LLf  •  •  •  / WN(1,2, '  •  •  '^^  d^  ' ' '   d3S  > 


\W  Z- 


*-uW2>)  r ...  rH"g(31;;';;g)  ^  ^  -  A 


\   ZN  ^   V  ^  (N-2)!\ 


2    r42)(£(g)^),3,  r      f»»#-'s)AA      a 


^    •••     ." J.   ^,_^    a^a-'a-  *jh 


-J-  ffn'2)(l(2'M)AA 


IV        *>rf 


X     ZN   V 


x  J  " " '  J  ,    ,  n   3(n-M     d  £  d  S  •  •  ■  d  2  + 


4=1  *   ~~  ZN 


Defining  the  activity  z  by 

dinZ 

M-\3z)  -  -  -^  ,  (2.3-7) 


dN   ' 


and  assuming  N  »  £   one  obtains 


in(\3z/  =  JnZ^  -  inZN  , 


Ik 


or 


J  -   X-3J  \L   .  (2.3.8) 


Thus,  taking  the  limit  as  N  -»  <»,  the  expansion  of  the  pair  distribution 
function  in  powers  of  the  activity  becomes 


n  (1,2)  =  z  Z  <«b(.2)(l,2)zi  .  (2.3.9) 


Using  Eq.  (2.3«*0  with  h=l  the  density,  n,  can  be  found  as  a 
power  series  in  z  also 


00 

n  =  E  SbQ7~   ,  (2.3.10) 

£=1     l 

where  we  drop  the  superscript  on  b,  (l)  and  omit  the  argument  since  we 
are  dealing  with  a  homogeneous  fluid  in  which  X>\       is  independent  of 
position,.   It  should  be  noted  that  homogeneity  also  requires  that  the 
2 -particle  cluster  integrals  depend  only  on  r   „    From  Eqs .  (2.3.9) 
and  (2. 3 '10)  a  series  for  n  (r, _)  in  powers  of  n  can  be  constructed 


it—  cL 


The  first  two  coefficients  of  the  series  are  given  by 
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y   (r  )  =  b^(r  ) 
r2K    12'    1   v  12y 


73(r12)  =  2i2\r12)    -   ^b(2)(r12) 


(2=3.12) 


The  function  obtained  experimentally,  usually  by  x  ray  or 
neutron  diffraction,  is  the  Fourier  transform  of  the  pair  correlation 
function 


,    ,    n2(ri2} 
S(r12)  =  g 


(2.3.13) 


n 


From  the  comments  following  Eqs .  (2.1. 3)  and  (2.1.4)  it  is  seen  that 
the  pair  correlation  function,  g(r   ),  is  normalized  so  that 


rli5ooS(ri2)  =  1   > 
12 


(2.3.1^) 


provided  that  the  interparticle  potential  has  a  finite  range.   If  one 
defines  &±(r12) >    i=l,2,3,...,  by 


:(r   )  =  b^(r   ) 
A  12'     1  K    12' 


1  +  Z  n1g  (r   ) 
i=l 


(2.3.15) 


it  is  seen  from  Eq.  (2.3.12)  that  the  coefficient  of  the  term  linear  in 
the  density  is 


5i(ri2}  -  A2),    ; 

bl      (rl2} 


2b^(r      )    -   kb  b^hr      ) 
2       ■    12y  2   1      ^    12' 


(2.3.16) 


16 


which  can  be  written  in  terras  of  the  diagonal  density  matrix  elements, 
W,  as 


gl(ri2}  "  W0(l,2)J 


.W3(l,2,2)  -  W2(1,2)W2(1,^) 


W_(l,2)W  (2,^)  +  W_(l,2) 


d3£.   (2.3.17) 


The  linear  term  in  the  density  expansion  of  g(r   )  (or  n  (r   ))  is  thus 
expressed  in  terms  of  the  density  matrix  elements  of  systems  of  three 
or  fewer  particles. 

Similar  manipulations  with  the  cluster  expansion  of  W  yield 
the  well  known  virial  expansion  for  the  pressure  P  (for  example  see 
Huang  [10]) 


Z^  _  !  +  B(Tl   C(Tl  (2  3  18) 


where  the  third  virial  coefficient,  C(T),  is  given  in  terms  of  the 
cluster  integrals  by 


=  -N?(2b_  -  4bQ2)  f  (2.3.19) 


'0V"3 


and  where  N  is  Avogadro's  number  and  V  is  the  molar  volume. 

Using  the  definition  of  the  cluster  integrals  and  of  g  (r  p),  C  can  be 

written  as 


IT 


N? 


f^12)^12^\-   kh22]      '  (2.3^29) 


Since  W  (r   )  and  b  have  already  been  evaluated  for  the  Lennard- Jones 
(12-6)  potential  [6]  it  remains  to  calculate  g  (r   )  in  order  to  obtain 
both  the  pair  distribution  function  to  first  order  in  the  density  and 
the  third  virial  coefficient. 
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3*    PATH  INTEGRAL  FORMULATION  FOR  THE  DENSITY  MATRIX 

3 .1   Expression  of  Density  Matrix  Elements  in  Terms  of  Wiener  Integrals 

The  expression  of  the  quantum  mechanical  density  matrix  in 

terms  of  functional  integrals  has  been  treated  by  several  authors 

[11,12,13,1^,15] •   ^ne  mathematically  rigorous  form  of  the  functional 

integral  in  question  is  the  Wiener  integral  or  expectation  value, 

E{F[x(t)]),  of  a  functional,  F[x'xt)],  over  Wiener  measure  on  the  set  of 

continuous  functions,  x(t),  over  a  finite  interval  in  the  independent 

variable  To  The  following  treatment  will  make  use  of  the  Wiener 

integral  formulation  [11,12,13,1^],  which  exhibits  normalization 

constants  more  explicitly,  rather  than  the  more  intuitive  notation  of 

Feynman  [15]>  which  clarifies  the  transformation  properties  of  the  path 

integral . 

The  quantity  directly  related  to  the  density  matrix  is  the 

conditional  Wiener  integral,  which  involves  Wiener  measure  on  the  subset 

of  continuous  paths  which  have  fixed  values  at  both  endpoints  of  the  T 

interval.   Let  the  3N-dimensional  vector  (the  number  of  dimensions  is 

arbitrary,  but  3N  is  specified  for  later  application  to  the  system  of  N 

particles  in  3  dimensions),  r(r)  =  (r,  !,t),  . . .  ,r  (t)  ),  be  a  continuous 

vector  function  of  the  variable  T  on  the  interval  0  <  t  <  (3.   Let  the 

interval  [0,(3]  be  divided  into  k  equal  intervals  of  length  At  =  p/k, 

and  take  the  division  points  to  be  0  =  T_  <  T,  <  .  •  <•  <  T,  .  <  T.  =  (3. 

0    1         k-1    k 

Let  r(T;k)  be  a  continuous  vector  function  of  t  consisting  of  linear 

segments  over  the  intervals  (t„  . ,T. ),    i  =  l,2,.^.,k,  with 

r.  =  r(x.;k),  and  satisfying  the  endpoint  conditions  r  =0,  r  =  R. 
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If  the  functional  F[t(t)]  is  defined  on  all  continuous  functions,  r(i:), 
satisfying  the  above  endpoint  conditions,  then  the  conditional  Wiener 
integral  can  be  defined  by  the  limiting  process 


00         00 

E{F[£(T)]|r(p)  =  R}  =  lira  /  ...  /  F[r(T;k)  ]d^, 


■00      -00 


(3.1.1) 


where  the  measure  d|~u  is  given  "by 

K. 


3N      2 
2 


-3Nk 


duR  =  (2it)  '"  exp(— )(2itAt)     exp 


k  (r.  -  r.  ,  Y 


1=1 


2AT 


n  d3Nr.  .  (3.1.2) 
i=l   ~X 


Note  that  the  measure  dun  is  normalized  so  that 

k 


duk  =  1  , 


(3-1.3) 


and  thus 


E{l|r(p)  =  R}  =  1 


(3.1**0 


The  relationship  between  the  conditional  Wiener  integral  and 
the  density  matrix  is  given  by  a  theorem  due  to  Kac  [13, lM"   It  states 
that  if  the  scalar  field  U(r)  has  a  finite  number  of  discontinuities 
and  is  bounded  below,  then 
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■BE 


-3N 


Z.Y.  (R)¥.(R')e   J  =  (2*p)    exp 


(B'-fi)' 
2P 


X  E  iexp 


B 


U(r(x)  +  R)dx 


r(B)  =  R'-Rf  , 


(3.1-5) 


where  Y.(R)  are  the  normalized  eigenfunctions  of 
J 


KsiVs>  + 


Ej-U(R)  JSTjCr)   =  0 


(3.1.6) 


If  one  transforms  the  variable,  R,  in  order  to  make  Eq.  (3° 1.6)  into  the 
time  independent  Schrodinger  equation  for  a  system  with  the  Hamiltonian 


ii   2 
H  =-E.  7T-  ST.    +   U(R), 
1  2m  ~i    y~' ' 


(3.1.7) 


and  makes  the  further  transformation  t  -»  t/|3,  £(t)  "~*  £(t)/nP,  which 
normalizes  the  T  interval  to  [0,1]  while  leaving  the  Wiener  integral 
invariant,  then  the  general  density  matrix  element  can  be  written 


<F|e     |R'>  =  exp 


-^(R'-R)2 


X  EX  exp 


•eiufej(T)+5r 


r(l)  =  "P  (R'-5) 


(3.1.8) 
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This  relationship  holds  only  for  the  Boltzmann  matrix  elements  since  the 
sum  in  Kac '  theorem  is  over  all  eigenstates,  regardless  of  symmetry* 

As  stated  in  the  references  [13,1^-],  Kac1  theorem  also  requires 
that  U(r)  be  such  that  the  eigenvalue  problem,  Eq.  (3-1.6),  yields  a 
discrete  spectrum.   If  the  system  is  enclosed  on  a  box  of  volume  V,  then 
U(r)  should  have  a  component  due  to  the  walls  of  the  box.   This 
component  certainly  causes  the  spectrum  to  be  discrete  and  gives  a 
natural  normalization  for  the  eigenf unctions .   However,  if  one  restricts 
the  temperature  to  values  larger  than  those  at  which  long  range  order 
phenomena,  such  as  Bose  condensation,  become  important,  then  the 
potential  of  the  box  contributes  nothing  to  the  problem  in  the  limit 
V  -*  oo.   One  may  neglect  this  component  of  U(r)  provided  one  considers 
the  wave  functions  of  the  continuous  spectrum  to  have  delta  function 
normalization  and  treats  the  sum  over  states  Z.  as  representing  a  sum 
over  the  discrete  states  plus  an  integral  over  the  continuous  spectrum, 

A  slight  reformulation  of  Eq.  (3. 1-8)  will  be  useful  in  the 
following  work.  Wiener  measure  on  r(f)  gives  rise  to  the  interpolation 
formula  [l6] 


r(T')(T"-T)  +  r(T")(T-T')       r   ,    ,  w  ,,   v-|l 
,(T)  =  -,_, +|  L  (T(^»T)  T>  j  2    ,   (3.1.9) 


where  t'  <  t  <  t"  and  where  £  is  a  3N-dimensional  random  variable  whose 
coordinate  random  variables  are  independent  and  normally  distributed 
with  mean  0  and  variance  1=   From  this  formula  it  can  be  seen  that  if 
r'(T)  obeys  the  endpoint  conditions,  r'(0)  =  r'(l)  =  0  then  both  r(x) 
and  r'd;)  +  : — (R'-R)  are  distributed  as  the  random  variable 

~  A   ~   ~ 
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^-f1  (5'-£)  +1[t(i-t)]2 


Equation  (3. 1.8)  can  then  be  written 


,  "PHNi 
<R  e     R'>  =  exp 


-4(5'  -g)2 


X  E-^exp 


-PJU 

0 


Wiit 


r(x)  +  R  +  t(R'-R)  dT 


r(l)=Oh  , 


(3.1.10) 


where  the  prime  on  r'(x)  has  been  dropped. 

A  useful  intuitive  picture  of  the  Wiener  integral  is  obtained 
by  thinking  of  an  imaginary  motion  of  a  system  of  N  particles  over  a 
unit  "time'1  interval,  0  <  t  <  1.  At  time  T  =  0  the  coordinates  of  the 
particles  are  given  by  R  and  at  1  =   1  by  R' .  The  quantity  (R'-R)   is 
thus  the  sum  of  the  squares  of  the  distances  between  the  initial  and 
final  positions  of  each  particle,  and  the  Wiener  integral  in  Eq.  (3.1.10) 
is  the  average,  over  all  possible  paths  of  the  motion,  of  the  Boltzmann 

_Q  [T 

factor,  e   ,  with  the  potential  energy  U  replaced  by  the  average  of  the 
potential  over  the  path.   The  paths  are  weighted  by  Wiener  measure  which 
is  related  to  the  average  over  the  path  of  the  total  kinetic  energy  of 
the  system. 

It  should  be  noted  that 


nS 


:(t) 


2* 
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is  a  measure  of  the  deviation  of  the  motion  from  the  shortest  straight 
line  path  from  R  to  R1,  and  as  the  classical  limit  is  obtained  by  letting 
\  ->  0,    the   deviation  vanishes  in  this  limit.  The  off-diagonal  matrix 
elements  vanish  in  the  classical  limit  since 


lim  exp 
\-*0 


-4(5'  -g)2 


ee  6(R'-R)  ,  (3-l.ll) 


where  5  is  the  Dirac  delta  function,  and  the  diagonal  matrix  element 
approaches  the  Boltzmann  factor 


-PH^       -PU(R) 
lim  <R|e     JR>  =  e    ~   .  (3-1.12) 

\-0 


3-2   Upper  Bounds  on  Off -Diagonal  Matrix  Elements 

Both  diagonal  and  off- diagonal  matrix  elements  for  Boltzmann 
statistics  are  involved  in  the  calculation  of  g  (r   )  for  a  Bose  or  Fermi 
system.   For  example,  in  the  evaluation  of  the  3_particle  term, 


-f3H 
<!>2,2|e    3|1,2,3>+, 


two  types  of  off -diagonal  Boltzmann  matrix  elements  arise.   One  involves 
the  exchange  of  coordinates  of  two  of  the  three  particles  and  is  of  the 
form 


2k 


■PH 


<^2,^h   3|g,l^>  =  exp   -4 


X 


(2"1)2  +  (1"2)2 


x  E^exp 


-0  /  U,  -^r  r(x)  +  (2,1,5)  +  T[(i*2,^)-(2,1,2)]  dT 


3  US 


r(l)=0 


(3-2.1) 


The  other  type  involves  a  cyclic  exchange  of  the  coordinates  of  all  three 
particles  and  is  typified  by  the  matrix  element 


-PH 


<i>2,3|e   3|2,2,1>  =  exp  -±A   (2-l)2  +  (^-g)2  +  (1-3) 


X 


EJexp  •pJU3f-^r(T)  +  (2,3,1)  +  T[  (1,2, ^-(2,3,1  )]  jdT 


r(l)=o|  .     (3-2.2) 


The  form  of  the  factors  multiplying  the  Wiener  integrals  in 
these  expressions  yields  a  significant  estimate  of  the  size  and 
temperature  dependence  of  these  off-diagonal  elements  for  intermolecular 
potentials  for  which  the  particles  are  characterized  by  strong  repulsive 
cores.   Suppose  the  3_body  potential  satisfies 

U3^~l'~2'~3^  =  °°'  if  rij  <  °     fQr  any  pair  i^ssl»2>^'  (3.2.3) 

The  Wiener  integrals  in  Eqs .  (3-2.1)  and  ( 3-2.2)  must  vanish  if  the 
initial  distance  between  any  pair  of  particles  is  less  than  a.   Hence 
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wherever  the  matrix  element  is  nonzero  the  factor  multiplying  the  Wiener 

P  P 
integral  is  smaller  than  exp(-2ira  /\   )  for  the  2-particle  exchange  term 

p  P 

and  smaller  than  exp(-3:rtcj  /\    )  for  the  3-particle  exchange  term. 

k 

For  He  ,  approximating  the  repulsive  core  diameter  by  the 

-ft 
parameter  a  =  2.56  X  10   cm  of  the  Lennard-Jones  (12-6)  potential,  one 

finds  that  the  larger  of  the  two  multiplying  factors,  that  for  the 

-•6T  o 

2-particle  exchange,  is  approximately  e     »   Thus  at  10  K  one  expects 

the  off -diagonal  elements  to  be  smaller  by  a  factor  of  at  least  .0025 

than  the  diagonal  element,  for  which  the  multiplying  factor  is  unity. 

-f3H2 
The  argument  for  the  2-particle  exchange  term  in  <L,2|e     !i>2>  is 

the  same  as  the  argument  for  the  similar  3-particle  term  in  Eq.  (3»2.l) 

and  yields  the  same  result.   Due  to  restrictions  of  computing  time  the 

numerical  scheme  to  be  developed  below  for  the  calculation  of  the 

diagonal  matrix  element  is  limited  in  accuracy  at  low  temperatures .   In 

all  cases,  the  above  bound  on  the  exchange  terms  shows  them  to  be 

negligible  compared  to  the  numerical  uncertainty  of  the  diagonal  term. 

Therefore,  in  subsequent  sections  the  exchange  terms  will  be  discarded 

and  the  diagonal  Boltzmann  matrix  element  will  be  used  as  if  it  were  the 

same  as  the  diagonal  element  for  a  Bose  or  Fermi  system,. 

The  above  argument  was  first  used  in  connection  with  the 

2-particle  exchange  term  in  the  second  virial  coefficient  [17].   The 

suppression  of  this  2-particle  exchange  term  with  increasing  temperature 

is  actually  somewhat  stronger  than  is  indicated  above.   This  is  verified 
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both  by  numerical  calculations  [5,6]  and  by  upper  and  lower  bounds 
obtained  analytically  by  Lieb  [18],  who  shows  that  at  high  temperatures 
the  term  is  asymptotically  proportional  to 


-2 
Qualitatively  the  increase  in  the  constant  multiplying  -X,   in  the 

exponent  comes  about  because  the  Wiener  integral  in  Eq.  (3-2.1),  as 

well  as  the  multiplying  factor,  decreases  exponentially  with  increasing 

temperature.   This  is  due  to  the  fact  that  the  set  of  paths  for  which 

the  "motion''  of  the  system  from  time  t  =  0  to  time  T  =  1  brings  the  two 

exchanged  particles  no  closer  than  a   from  each  other  for  any  T  is 

assigned  exponentially  smaller  Wiener  measure  as  T  -»  <». 


3-3   Approximation  of  the  Diagonal  Matrix  Element  by  a  Multiple 
Piemann  Integral 

The  calculation  of  the  diagonal  matrix  elements  for  Boltzmann 

statistics  is  now  reduced  to  the  evaluation  of  a  Wiener  integral.   For 

the  diagonal  element  of  the  3-particle  density  matrix  Eq.  (3*1-8) 

becomes 
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<l,2,3|e   *\l,2,$> 


1 


=  E  -s  exp 


-Pj  U3  h=r(T)  +  (l,2,3)jdi 


r(l)=0  f  .  (3-3.1) 


Analytical  evaluation  of  the  above  Wiener  integral  is  possible,  however, 
only  in  the  cases 


U3(rx,r2,r  )  =  0  , 


and 


3 
U3(r1,r2,r3)  =  .^^r.)  , 


where  f.(r.)  is  quadratic  in  each  coordinate  of  r. .   Thus  for  a  realistic 
1  ~i'  ~i 

interparticle  potential  it  is  necessary  to  rely  on  a  numerical  evaluation 
of  the  Wiener  integral. 

Numerical  schemes  for  evaluating  Wiener  integrals  have  been 
proposed  and  applied  by  several  authors  [6,19,20,21,22,23]*   The 
simplest  of  these  methods  involves  taking  a  finite  but  large  value  of  k 
in  the  limit  definition  for  the  conditional  Wiener  integral,  Eqs .  (3.I.I) 
and  (3-1.2),  and  using  the  result  of  a  numerical  evaluation  of  the 
k-fold  multiple  integral  as  an  approximation  to  the  Wiener  integral. 
Since  k  is  large  and  the  weight  function  du  can  be  conveniently 

K. 

expressed  in  terms  of  Gaussian  density  functions  the  evaluation  of  the 
k-fold  integral  is  usually  done  by  the  Monte  Carlo  method. 
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This  simple  scheme  is  adequate  for  one  or  two  particles  under 
the  influence  of  a  computationally  simple  potential  and  has  been  applied 
successfully  to  such  problems  [6,20,22],     As  the  evaluation  of  the 
integrand  of  the  multiple  integral  becomes  more  complex,  however,  the 
computing  time  requirements  become  unreasonably  large,  and  for  this 
reason  methods  have  been  proposed  for  reducing  the  multiplicity  k  of  the 
integral  to  be  performed  [19,21,23] .   The  present  calculation  makes  use 
of  a  slight  modification  of  the  procedure  developed  by  Fosdick  [21]. 

The  Wiener  integral  in  Eq.  (3«3«l)  can  be  written  as  a  Riemann 
integral  of  a  product  of  Wiener  integrals  involving  a  smaller  value  of  3 
(higher  temperature)-   If  the  interval  0  <  t  <  1  is  divided  into  k 
equal  subintervals,  (t.  -,,t.),  i=l,2,...,k,  then  the  integrand  in 
Eq,  (3.3=1)  can  be  written  as  a  product  of  functionals  depending  on 
r(x)  over  only  one  subinterval 


F[r(i)]  =  exp 


f  u3  (^  £(t)  ♦  (i,s,3)j  **  ] 


n  fJrfT)]  , 
i=l 


(3-3.2) 


where 


fiteC-O] 


=  exp 


T. 

1 


-PJ   UJ_\_r(T)  +  (1,2, 3)  dT 


T.  ,   V2Jt 
l-l 


(3.3.3) 
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If  each  of  the  k  subdivisions  is  now  further  subdivided  into  n  equal 
parts,  Eq.  (3-1-1)  may  be  written  in  the  form 


00  00 


s 


{F[r(T)]|r(l)   =  0}    =  lim      /  ...    /  dunkF[r(T;nk) ]    , 

n->oo  ^  ^ 


•00  -00 


00  00 


00  00 


du        IT     lim 
i=l  n-»oo' 


/•••J\,nfi[^nk>]' 


00       -  CO 


-co     -00 


(3.3.10 


where  d(i  is  defined  as  in  Eq.  (3-1.2)  in  terms  of  the  coarse 

K. 

subdivision  only  and 


du,  „  =  (24)     2 


i,n 


exp 


(r.   -r.  _)' 
2(p/nk) 


3Nn 


X  exp 


■  n   (r.  .-r.  .  , Y 

-y     '~i;j  ~i;.]-iy 
.j=1     2(p/nk) 


(3.3-5) 


where  r.  .  =  r(f  .  .  :nk).   The  limit  on  the  right  hand  side  of 
Eq.  (3-3-M  is  just  the  definition  of  a  conditional  Wiener  integral 


00     00 


lim  /  .  . 
n-»  00^00 


du^f^r^nk)]  =  E^exp 


-p     U, 


Ti-1 


\ 


r(T)  +  (1,2,3)  Ut 


r(T.  .)  =  r.  ,,r(T.)  =  r.X   .     (3-3-6) 

~    1-1       ~X-1'~N  X  ~1J 
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The  time  interval  in  Eq<  (3«3-6)  may  be  renormalized  and  the 
endpoint  conditions  altered  by  the  transformations 


t1  =  k 


(x  -  x.^)  , 


'  (t)  =vk[r(i)  -r.  ,  -  i(r.  -r.  n)l 
~      -l-l    x~i   ~i-l/J 


Equation  (3«3-^)  then  becomes 


00         00 


E(F[r(T)]|r(l)  =  0}  = 


dn, 


(3- 3- 7) 


-00      -  X) 


X  H  Ei  exp 
i=l 


.& 


r(x) 


U-,  -=r  ~T=~  +  r-  1+T(r.-r.  J 


+  (i*2,3)WT 


r(l)=0^  , 


(3.3-8) 


where  the  primes  on  r(x)  and  t.  have  been  dropped.   Since 

\//k  =-  2rth2(p/k)/m  2  , 


the  Wiener  integrals  on  the  right  of  Eq.  ( .3 •  3 •  8 )  involve  a  temperature 
k  times  that  in  the  Wiener  integral  on  the  left.  We  proceed,  then,  by 
finding  a  high  temperature  approximation  for  the  Wiener  integral  and 
using  Eq.  (3«3'8)  "to  extend  its  range  of  validity  to  lower  temperatures. 
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The  high  temperature  approximation  to  be  used  is  the  well 
known  Wigner-Kirkwood  expansion  in  powers  of  |3  and  \  [1,2^].   The  most 
useful  form  of  this  expansion  for  the  following  work  will  be  the 
expansion  for  the  quantity 


G  = 


inEjexp  -p  /  U3  I  —  £(t)  +  r  +  T(s-r)jdi 


■(l)=0 


(3.3-9) 


Introduce  a  parameter  t]  multiplying  the  path,  r(l) ,    in  the  Wiener 
integral,  and  define  the  following  functions  of  r\ 


f(ti)  =  exp 


P  /  Un  M=L  r(T)  +  r  +  T(s-r)]dT 
J      3  V^  ~      ~     ~  ~  J      ■ 


F(tj)  =  E(f(ii)|r(l)  =  0}  , 


and 


G(t])  -  ZhF(r\)    > 


(3.3.10) 


A  formal  series  expansion  for  G  may  then  be  found  by  expanding  G(t|)  as 
a  power  series  in  tj  about  t)  =  0  and  taking  the  limit  r)  -»  1  to  obtain 


G  =  lim  G(t))  . 
1^1 


(3-3.11) 
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Up  to  order  t)  one  finds 


G(n)  =inF(0)  +^T^ 


F(0)   drj 


q=0 


2 


1   dF(n) 


_F(0)   dr] 


"2    1  d2F(n) 

n  J    F  0)   ,  2 
T]=0  -1     '    dr) 


T)  =  0 


(3-3.12) 


The  integrand  of  the  Wiener  integral  F(0)  is  independent  of  the  path  so 
it  is  seen  from  the  normalization  condition,  Eq.  (3.1.U),  that 


F(0)  -  f(0)  =  exp 


■P  /  U  (r  +  T(s-r))di 


(3.3.13) 


The  first  term  in  the  expansion  for  G  is  then 


G  -  -0  /  U  (r+T(s-r))dT  , 


0 


(3.3.1M 


which  is  just  the  exponent  of  the  classical  Boltzmann  factor  with  the 
potential  replaced  by  its  average  over  the  straight  line  path  from 
r  to  So   If  this  approximation  is  used  on  the  right  hand  side  of 
Eq.  (3- 3»8),  the  result  is  precisely  that  of  using  Eq.  (3.1.1)  with  a 
finite  value  of  k  as  mentioned  above. 
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The  coefficient   of  r\   in  the  expansion  for  G  is 


1       dF(n) 
F(0)     d-q 


n=o 


F(0) 


E^-pf(O) 


0 


-y —  r(i )    •    W   (r+T(s-r))dT 
>/2rt   ~  ~  3  ~ 


r(l)=o|  ,  (3-3.15) 


where  V  is  the  9-dimensional  gradient  operator.   If  one  reverses  the 
order  of  the  Wiener  integration  and  the  integration  over  1,    the  Wiener 
integral  depends  on  r(x)  only  for  a  fixed  1.      This  Wiener  integral  may 
be  evaluated  by  using  the  interpolation  formula,  Eq.  (3 .1.9)*  f°r  r('T') 
with  t'  =  0,  t"  =  1  and  taking  the  expectation  value  over  the  normal 
distribution  of  the  random  variable  £ 


E{r(T)|r(l)  =  0)  =  <|[t(1-t)]   >   =  0  , 

~  d  V 


(3.3.16) 


therefore 


dF(r,) 


4l 


=  0 


TJ=0 


(3.3.17) 


The  coefficient  of  r\   thus  vanishes  and  the  coefficient  of 


t)  1 2.   becomes 


3^ 


1  _2 


dn 


n=o 


r(l)=0 


2  9  9 


o- 


+  TP  J  ST  *  -Z,  ri(T)rJ(T)  5I5J  U3(^T(^r)) 


dT 


i=lj=l 


r(l)=Oh  , 


(3-3.18) 


au. 


where  r.(i)  is  the  i-th  coordinate  function  of  r(x)  and  vr^  denotes 
partial  differentiation  of  U  with  respect  to  its  i-th  argument.   The 
first  term  on  the  right  of  Eq.  (3.3.18)  is  of  order  0  \  «T    while 
the  second  term  is  of  order  (3\  «T   .  Thus  to  lowest  order  in  l/T 


1   d^F(n) 


T]=0 


9   9   .2 


~ "  /   ^   ^  ^|U  U  (r+T(s-r))E{r.(T)r.(T)|r(l)=0) 
1=1  .1=1    J  3  J 


2* 


dT  . 


0 


(3.3.19) 


Using  the  interpolation  formula  again  one  finds 


Efr, (T)r  (f)|r(l)  =  0}  = 


0  ,       (ifcj) 
t(1-t)  ,   (i=j) 


(3.3.20) 


Therefore  to  order  T 
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G  - 


r  2  r 

-P  /  U  (r+T(r-s))dx-^-  /  V  U  (m(s-r)h(l-T)di 


(3-3.21) 


3 


If  Eq.  (3.3-21)  is  used  for  the  Wiener  integrals  on  the  right 
of  Eq.  (3-3-8)  one  obtains  the  approximation 


00     00 


K{F[r(T)]|r(l)=Oj  =  J  .  . .  J  dnkF[r(-r;k)  ]Ck  ,        (3-3-22) 


-00    -00 


where 


C  =  n  exp 
k   1=1 


^i-l^^i^i-l) 


+  (1,2,3)1  T(l-T)dT 

(3.3.23) 


If,  for  each  i=l,2,...,k,  the  exponential  in  Eq.  (3-3-23)  is  replaced 

by  the  first  two  terms  in  its  power  series  expansion,  it  is  seen  that 
-2 


to  order  k 


Ck  =  n 
k   1=1 


1 


B\~ 


kitti 


F 


IL 


\ 


Hn& 


r.  , +f(r . -r .  . ) 
~i-l   v~i  ~i-l' 


+  (1,2,3)  T(l-r)dT 


(3-3- 2k) 


Equation  (3.3.22)  with  Eq,  (3-3.24)  for  C   is  precisely  the 

formula  obtained  by  Fosdick  [21]  who  shows  that  if  U  is  sufficiently 

-2 
regular,  the  formula  is  correct  to  order  k   .   In  the  present 

3 
application,  where  U  (1,2,3)  =  £  V(r..)  with  V(r)  given  by  Eq,  (2.1.2) 

the  regularity  conditions  are  not  satisfied  due  to  the  singularity  of 
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V(r)  at  r  =  0,   It  has  been  found  from  calculations  of  the  2-particle 
density  matrix  elements  that  the  use  of  Eq.  (3-3.23)  -^or  C,  S^ves  better 

K. 

results  for  small  k  than  are  obtained  using  Eq.  (3- 3*24)  due  to  the  fact 

that  at  the  singularities  of  U  ,  Eq.  (3»3«23)  is  bounded  in  absolute 

value  while  Eq*  (3« 3 -2k)    is  not.   For  this  reason  the  exponential  form 

for  C,  will  be  used  in  the  following  calculations, 
k 

The  Wigner-Kirkwood  expansion  can  be  carried  out  in  the  above 
manner  by  repeated  use  of  the  interpolation  formula  to  obtain  approxi- 
mations which  are  correct  to  order  k   for  any  m  under  the  appropriate 
regularity  conditions  on  U.   However,  the  increasing  complexity  of  the 
C,  obtained  in  this  way  causes  an  increase  in  computing  time  which  must 
be  balanced  against  the  decrease  obtained  by  using  a  smaller  value  of 
k  in  Eq.  (3.3*22).   The  special  case  k  =  1,  which  is  equivalent  to 
expanding  the  logarithm  of  the  diagonal  matrix  element  directly,  is 
simplified  by  the  fact  that  the  1   integrals  may  be  done  immediately. 
The  result  of  the  expansion  to  order  T   in  this  case  is 


in  Eiexp 


-^U3fj^r[r(T)+(l,2,3)AdlT  r(l)=0J 


2  2.  2 

-pu3(i,2,a)  -§^^3(1,2,3)+  |£-[vu3 (1,2,2)]' 


-^vV  (1,2,^)  '  (3.3.25) 

96O7T      ^ 
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k.        EVALUATION  OF  g^r   ) 

l+.l   The  Multiple  Integral  Approximation  for  g  (r   ) 

The  complexity  of  the  3_Particle  functions  entering  into  the 
calculation  of  g-,(r1p)  can  be  reduced  somewhat  by  consideration  of  the 
symmetries  of  the  problem.   Since  there  is  no  fixed  external  field,  the 
density  matrix  must  depend  only  on  the  relative  positions  of  the 
particles,  and  hence  the  first  simplification  should  be  to  separate  out 
the  center  of  mass  motion.   Define  the  usual  center  of  mass  and  relative 
coordinates  for  a  system  of  three  particles,  Figure  1,  by 


£CM  =  fvi-2-1)  > 


E^-g-A  . 


S  =  3-|{g+l)  .  (k.l.l) 


The  Jacobian  of  this  transformation  is  unity,  and  furthermore 


£j^-£(!*.*«£*iii  <*•"> 


so  that  the  coordinates  r_„  r,_,  and  Q  behave  as  independent  particles 

~CM  ~12      ~ 

of  masses  3^,  m/2,  and  2m/ 3  respectively. 

Taking  U  to  be  the  sum  of  pair  interactions  as  in  Eq.  (2,1.1) 

1  1 

and  noting  that  r,  _  =  3-1  =  Q  +  o  rio  an(3  r_„  =  3~2  =  Q-—  r,_  it  is  seen 

~13    ^  ~   ~  d       -L^-  ~2J    ~  ~   ~  2  ~12 
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FIGURE  1.   CENTER-OF-MASS   AND  RELATIVE 
COORDINATES  FOR  3  PARTICLES. 
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that  W  (1,2,3)  can  ^e  separated  into  the  product  of  a  center  of  mass 
factor  which  is  independent  of  U  and  a  factor  containing  only  the 
relative  coordinates,  r,   and  Q, 

W_(l,2,£)  =  <l,2,2|e  '  *|l,g,3>  , 


"PKPM.  ,    "PHrpl, 

=  <We  l£cM><£l2^h  |r12,§>,  (4.1.3) 


where 

K  -  *% 

CM    2 (3m) 
is  the  kinetic  energy  of  the  center  or*  mass  motion  and 

Hrel  "  -2^27  "  iWT)  +  U3(^12'S)  &**) 

is  the  Hamiltonian  operator  for  the  relative  coordinates.   It  is  to  be 
understood  that  the  normalization  factor  of  X   for  each  particle 
coordinate  appearing  in  the  definition,  Eq_.  (2.2.4),  of  the  Boltzmann 
matrix  elements  is  to  be  replaced  by  the  thermal  wave  length  for  a 
particle  of  mass  corresponding  to  that  associated  with  the  appropriate 
transformed  coordinate 


(4.1.5) 


ko 


The  normalization  of  the  density  matrix  elements  yields 


<j"    e 

-CM1 


•pKc: 


M 


r  >  =  1 
~CM 


(4.1.6) 


so  that,  using  Eq.  (3-1=8),  the  diagonal  density  matrix  element  can  he 
expressed  in  terms  of  a  Wiener  integral  over  paths  in  the  6-dimensional 
relative  coordinate  space  of  r, _  and  g: 


W3(i'2,2) 


EJexp  -Pj  U3(r(T),q(T);r12,3)dTj  r(l)=q(l)=o|  (^.1-7) 


where 


U  (r(T),a(T);r   ,g)  =  V  —  r(T)  +  r 


3 


^T 


-12, 


12 


+  V  -^-  (nT?  ^r)+r(r))+^+Z2-  )  +  V  -~  (^3  ^T)-r(x))+^ 


-12 


2^7 


2</jT 


(4.1.8) 


and  V(r)  is  given  by  Eq.  (2.1.2). 

The  center  of  mass  motion  can  also  he  eliminated  from  the 

2-particle  density  matrix  elements  occurring  in  the  expression  for 

g  (r.p),  Eq.  (2.3-1?)  •   A  typical  2-particle  diagonal  matrix  element  is 


in 


W_(l,2)  =  E^exp 


J=T3,(t)  +  S  +  ^ 


a'(i)=o 


C^.i.9) 


One  can  thus  write 


=  l(ri2)  =  OlTiy/  d3§ 


EJF123(r(T),a(T):.|,E(l)=a(l)=0 


-B|F12(r' (t)  )F13(q' (t)  )|r '(!)=£' (l)=oj 


-BJF12(r\'T))F23(q:;(T)) 


r"(l)=aM(l)=0}-+  E^F12(r'"(T))  r"*(l)=0 


(^.1.10) 


where 


F123(r(T),g;(T))  =  exp 


-P  /  U  (r(T),a(T);r   ,g)dT 


3 


1 


F12(r*(T))  =  exp 


\ 


0   '"* 


1 


F13(a'(T))  =  exp|-pj  Vlj^'CO  +  S  +  ^ldT 


and 


F23(q"(-))  =  expf-p  jVj^q"(T)  +  Q 


-12 


dT 


tt.l.ll) 


Since  the  expectation  value  is  a  linear  operator,  one  may  identify 
r'(T),r"(T),  and  r'"(T)  with  r(i)  and  g/ (  t  )  and  &"(t)  with  q(x).   The 


quantity  g  (r   )  can  then  be  written  in  terms  of  a  single  Wiener  integral 
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Sl(rl2}  =  W  (1,2)  /d3S  E{F123(~(T)^(T))"F12(£(T))F13(a^)) 


-F12(r(x))F23(q(T))+F12(r(T)) 


Now  let 


r(l)=q(l)=Ol  .  (4.1.12) 


G(r12,S)  =  w  (^2)  W3(l,2,3)-W2(l,2)W2(i,3)-W2(i,2)W2(2,3)4W2(l,g) 

(fc.1.13) 

be  the  g- integrand  of  g  (r   ).   Noting  that  the  interaction  potentials 
depend  only  upon  the  magnitude  of  the  separations,  r   ,  r-io>  anc*-  rpo 
between  the  particles  it  can  be  seen  that  if  r   is  taken  to  lie  along 
the  axis  of  polar  coordinates,  (pn,  6n,    cp  ),  for  Q,  then  G(r   ,Q)  is 
independent  of  the  direction  of  the  polar  axis  and  of  cp  .   Thus  one  can 


write 


Sl(r12)  =  SnjTsin  e^f^   C(r12,pQ,eQ)  .     (4.1. 14) 


The  integral  over  the  position  of  the  third  particle  which  appears  in 
g  (r   )  can  thus  be  reduced  to  a  twofold  integral. 


Equations  (4.1.8)  and  (4.1.11)  show  that  F    can  be  written 


as  a  product 


F123(r<T),a(T))  =  F12(r(T))F13  (^q(T)+|r(x)  F^  [^Sft(T)^.(T)|  . 


(4.1.15) 


^3 


It  should  be  noted  that  if  r(x)  and  £(t)  are  conditional  Wiener  paths 

1 

2 
with  r(0)=£(0)=r(l)=q(l)=0  and  variance  [t(1-t)]  ,  then  the  linear 

combinations 


2  %(?)  ±  ^(t) 


are  conditional  Wiener  paths  with  the  same  distributions  as  r(x)  and 
£(t).   This  can  be  seen  from  the  interpolation  formula,  Eq.  (3-1-9)-   If 
we  consider  just  one  component  of  the  3-dimensional  paths,  the 
interpolation  formula  gives 

1  1 

q  (t)  =  £  [t(1-t)]2  ,  r*(T)  =  i Jt(I-t)]2 


and 


^   qx(T)  ±   |rx(T)  =(&   6q  +  |  U[t(1-t)]2  ,        (4.1.16) 


where  £   and  £   are  normally  distributed  independent  random  variables 
q      r 

with  mean  0  and  variance  1.   The  linear  combinations 


JT      1 
—  £   +  —  £ 
2  5q  -  2  5r 


are  thus  normally  distributed  also  and  have  mean  0  and  variance 

3   1 

j-  +  j-  =  1.   The  distributions  for  the  coordinates  of  q(x),  r(i), 

and 
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^£(T)  ±  \  £(T) 


are  thus  equivalent. 

Define  now  the  quantity 


gc(r12)  = J   d3Q  Gc(r12,Q)  ,  (^.1.17) 


where 


Gc(r12>^  =  W^TTU  4F123(r(T)^(T))-F12(r(T))F13p  q(r)  +  \   r(T) 
"  F12(r(T))F23^£(T)-|r(T)lF12(r(T))jr(l)=3(l)=o} 


(4.1.18) 


Due  to  an  oversight  the  quantity  g  (r   )  rather  than  g(r   )  was 


originally  calculated.   The  difference  between  these  two  quantities 
arises  because,  although  the  distributions  of  the  arguments  of  F,   and 
FQ_  in  Eq.  (4.1.18)  are  the  same  as  the  distributions  of  the  arguments 
of  these  functions  in  Eq.  (4.1.12),  the  use  of  a  linear  combination  of 
jj(t)  and  r(x)  as  the  argument  of  F   in  Eq.  (4.1.18)  introduces  a 
spurious  correlation  between 

F13fl(')  -  |  E(t) 

and  the  functional  F   (r(r))  multiplying  it  and  similarly  for  the 
product  F,   F   o   The  two  factors  in  these  products  should  be 
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uncorrelated  since  they  arise  from  the  product  of  two  independent 
Wiener  integrals.  The  function  g  (r   )  can  be  written  as  the  sum  of 
g  (r   )  and  a  term  which  corrects  for  the  unwanted  correlation: 


where 


and 


gl(r12)  =  gc(r12)  +  ge(r12)  , 


*e(r12>  =/d3S  Getr12,§)    , 


(4.1.19) 


Ge(ri2^)  =  w^yE{Fi2^))k3fFa(^  +|£(t)   -  f13(s(t)) 


+  F23PT^)  -  |r(T)J  -  F23(q(T)) 


r(l)=3(l)=o| 


(4.1.20) 

The  numerical  calculation  will  thus  involve  separate  evaluations  of  the 

terms  gQ(r12)  and  §e(ri2^ 

Using  Eq.  (4.1.15)  the  quantity  G  (r      ,Q)   can  be  written  in 


cv12' 


the  form 


°c<ri2'S>  =  iOT^)  E{Fi2(£(l))[ri3lT  s(^) + 1  iM 


-   1 


f23Pf^t)-|^t)I-  1 


r(l)=q(l)=0 


(4.1.21) 


The  results  of  Section  3-3  niay  now  be  used  in  order  to  express 
G  (r,p,Q)  and  G  (r   ,Q)  in  terms  of  multiple  Riemann  integrals  over  a 
6k-dimensional  space.   Let  r(T;k)  be  the  broken  straight  line  path 
defined  in  Section  3-1  with  breakpoints  at  r.  /i=0,l, . .  ..,k,  and  similarly 
for  q(x;k).   Define  the  functions 


k6 


D12(r(T;k))  =  F12(r(T;k)) 


)   -  Z     -&£   I  V"[^;[r,    ,«(*,-&■•.,  )]    +r1Q|T(l-T)dT 
L  1=1  2nk 


V?"1- 


■i  ~i-l/J     ~12 


0 


D13(a(T;k))    =   F13(q(T;k)) 


1 


2  r, 

2*k  0 


£12 


exP^    &L.J    V"  ^[ai.1+T(a.-ai.1)]^+^  T(l-T)dT 


D23(£(T;k))   =  F23(q(x;k)) 


X  exp 


k     M2 


L^/vWa-+^-3-)]+s- 


-12 


T(l-T)dT 
J 


where 


V"(r)    =  V2  V(r) 


(V.1.22) 
(4.1.23) 


The  approximations   for  G     and  G     are  then 

c  e 


00  00 


Gc(r12,g) 


Vi'S) 


D12(r(T;k)) 


D13 \~2    ^T;k)  +  \  ~(T;k) 


-   1 


-00  -00 


X      D 


23 


Sa(T;k)   *      (T;k))-1 


rk     qk 


(1+.1.24) 


and 
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00         00 


Ge(r12'S)  ~  OT^y  /•••/D12(r(x;k)) 


-oo      -oo 


d13I  ^  q(x;k)  +  -  r(T|k) 


D13(q(T;k))  +  D23l  —  q(T^k)  -  ^  r(xjk)  I    D23(q(T;k))|  d^rkd^  k  > 


(4.1.25) 


'3 


where 


3     _3k 
dngk  =  (2«)2(2*/k)  '   exp 


f   k   (s  ,-s   J2"]^! 


(4.1.26) 


with  s  =  r  or  q, 

From  the  discussion  in  Section  3-3  of  the  order  of  accuracy 
of  Eq.  (3»3»22)  as  an  approximation  for  the  Wiener  integral  it  is  seen 

that  the  above  approximations  for  G  and  G  cannot  be  expected  to  be 

-2 
correct  to  more  than  order  k   .   The  order  of  accuracy  may  be  less  due 

to  the  singularity  of  V(r).   In  view  of  this  fact  the  integrals  over 

T  in  Eq.  (4.1.22)  may  be  replaced  by  their  Simpson's  rule  approximations 

without  decreasing  the  order  of  accuracy  of  these  formulae.   If  we  let 

r    and  q    be  the  3k-dimensional  vectors  (r, ,  r0,,°>,r,   )  and 

(qn  )    <i^)-'')(li   )  respectively  and  recall  that  r„  =  q_  =  r,  =  q,  -  0, 
~17  22.  ~k  ~0   ~0   ~k   ~k 

then 


k8 


00  00 


Gc(r12'S) 


w,<i'g) 


°12(£(k)> 


,s/T     (k)  1  .  (k)' 
D13      2     S       +2* 


-00  -00 


'] 


»J#»""-h""U 


d^rk%k   ' 


(U. 1.2.T) 


and 


00  00 


».('«.«)  ■  obr/  •  •  •  /  ^(k) >W?  *(k)  +  I  ;(k) 


-  Di3<a(k)) 


-00  -00 


♦"J?»W-*sW>)-%3<a(k>> 


du  .  d(j.  .     , 
rk     qk 


(i^.1.28) 


where 
D12(£(k))   =   < 


{ *  iBfe  ** 


£12  I  1^  v       2  ;        -12 


1  A,  I        \       _ _ ? t /  A,     /~i  ~i-l\ 


23J 


Hil.21 


\ 


-12' 


l^^i.i+3tTl^vl^( 


X      2i+Si-lN  „  £12 


)+Q+- 


i=lu  W« 


*2       X2     Jx  ,V*I-1 


+  VI  —  Q.+Q+  -~|  +  £—  V'!   —  (- 

wit 


0+Q+ 


-12 


2-  2 


(^.1.29) 


where  the  +  sign  in  the  last  equation  corresponds  to  D   and  the  -  sign 


to  D   . 


h9 


4.2   Monte  Carlo  Evaluation  of  the  Multiple  Integrals 

It  has  "been  indicated  above  that  the  6(k-l)-fold  integrals 
appearing  in  Eqs .  (4,1 ,27)  and  (4 „ 1,28)  can  be  evaluated  by  a  Monte 
Carlo  procedure o   The  quantities  d(i   and  du   have  the  properties 
of  probability  densities  on  3 (k-l) -dimensional  spaces  since  the 
integral  of  either  of  these  quantities  over  any  3 (k-l) -dimensional 
Borel  set  is  non-negative  and  the ' integral  over  the  whole  space  "is 
unity  (see  Eq<-  (3«1»2)).  The  simplest  Monte  Carlo  scheme  for  eval- 
uating the  multiple  integrals  is  to  choose  a  sequence  of  independent 

(k)   (k) 
6k-dimensional  vectors  (r   ,    q   ),  r    =  q    =  0,  i  =  1,  2,  . . . ,  M 

~i    ~l  ~k,i   ~k,i 

with  probability  density  du   du  ,  evaluate  the  integrands  of  G  and 

G  for  each  vector,  and  average  over  the  set  of  M  vectors-   If  one 

lets  I  (r  ■   ',  q   )  and  I  (r    ,    q   )  represent  the  integrands  of 

G  and  G  respectively,  then  by  the  law  of  large  numbers 


U»  i  Z     I(r(k)   q(k)} 


-/•■■■/  Wk)»3(k)) 


(4,2,1) 


dn  .  djj.   =  G  } 

rk   qk    c 


-00  -00 


with  probability  one  and  similarly  for  G  , 
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(k) 

Vectors  r    with  distribution  du  .  can  "be  generated  from 
~  rk 

3(k-l)  independent,  normally  distributed  random  variables  by  means  of 

the  interpolation  formula,  Eq.  (3«1.9) •  The  endpoint  conditions 

require  that  r  =  r.  =  0,  and  the  other  3(k-l)  coordinates  can  be 
^         ~o   ~k  v    ' 


generated  by  the  equation 


:i-i  (1  -  l> 


+    k 


-  (i  -  -) 

k  v    k' 


1 


i-1 


1 

-I  2 


(4.2.2) 


for  1  a  1,  2,  . ..,  k  -  1,  where  i.    is  a  3 -dimensional  random  variable, 

the  coordinate  variables  of  which  are  normally  distributed  with  mean  0 

(k) 
and  variance  1.   The  vectors  q    are  generated  by  the  same  procedure. 

~J 

This  straightforward  Monte  Carlo  scheme  has  been  used  [6]  to  compute 

the  2-particle  density  matrix  element,  Wp(l,2). 

The  3-particle  functions  G  and  G  are  well  adapted,  to  a  .  n 
more  sophisticated  Monte  Carlo  procedure  which  should  have  a  smaller 
variance  than  the  straightforward  scheme,  thereby  reducing  the  number 
of  samples  necessary  for  a  given  accuracy.   This  procedure  is  the 
importance  sampling  technique  first  introduced  by  Metropolis  and 
collaborators  [25]  for  evaluating  averages  over  the  Boltzmann 
distribution.  The  following  discussion  is,  to  some  extent,  an 


P-L 


extension  to  a  continuous  sample  space  of  the  description  of  this 
technique  given  by  Hammersly  and  Handscomb  [26]. 
Let  G  and  G  both  be  typified  by 

00         00 


where 


^'^'-v^r'^'s^         ^-h) 


and  let 


P(r<k»)  =  V^  (^-5) 

W,(l,2) 


This  factor,  vhich  appears  in  the  integrands  of  both  G  and  G  ,  is 
the  basis  for  the  importance  sampling  procedure.   It  is  seen  from 
the  definition,  Eqs  ~  (4,1,29),  of  D,„(r  ')  that  this  quantity  is 
non- negative*  Furthermore,  to  the  order  of  the  approximation  used 
in  Eqs,  (4<,1,24)  and  (4,1,25),  one  can  write 


00         00 


w2(i,2)  =/..-/  d^Ce00)  *nrk  >  e^.2.6) 
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and  hence 


00         00 


f,.,fl>(rM)   dnrk  -  1  .  (4.2.7) 


-00        -00 


(k) 
Thus  P(r   )  du   has  the  properties  of  a  probability  density  on  the 
~       rK 

3 (k-l) -dimensional  space  of  r   .   If  vectors  r    and  q^   are 

generated  with  probability  density  P(r    )  dp:   du   and  L(r   ,  q    ) 

^       ric   q^      '**'     ^ 

is  averaged  over  the  vectors  so  generated,  the  average  will  approach 
G  with  probability  one  as  the  number  of  samples  becomes  infinite. 
Qualitatively,  the  advantage  of  this  scheme  over  straightforward 
sampling  is  that  some  of  the  variation  in  the  integrand  l(£  /  q   ) 

has  been  absorbed  into  the  probability  density.   The  variance  reduction 

(k) 
is  indicated  intuitively  by  the  fact  that  vectors  r    for  which 

(k)  (k)   (k) 

D  p(r   )  and  hence  l(r   ,  q   )  are  small  will  be  generated  with  low 

probability  by  the  importance  sampling  procedure. 

(k) 

It  remains  to  show  how  vectors  r    can  be  generated  with 

(k) 
probability  P(r   )  du   .   This  is  accomplished  by  producing  a  Markov 
~      rK 

(k) 
chain,  {r   ,  t  =  1,  2,  3>  •  -  - } ;  with  a  stationary  transition 


•obability  function  p(r^  ,  |r|   ),  satisfying 
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J 


../p(r'W)d„ 


00         00 


•f-If-  JW(k)ii-<k)>  p(^(k>)  d^kd3(k_1)£'  -     ^-2-8> 


-00        -00 


where 


d3(k_l)r'  =  d3r'd3r'  ...  d3r"    ,  (^.2-9) 

~1  ~2       ~k-l 


and  where  J  is  any  set  in  the  3(k-l) -dimensional  space  of  r1    which 

is  measurable  with  respect  to  P(r*    )  du  ,  .   The  law  of  large  numbers 

~       r  k 

(k) 
for  Markov  chains  [27]  states  that  if  the  Markov  chain,  {r \    )>    is  such 

that  Eq.  (U.2.8)  is  satisfied  and  if  there  is  only  one  ergodic  set  then 


M  CO         00 

£a£^s(k)>-/---JVk). 


M~  M  ",  -«'"'  ~^"%'  f-f-**'   S^''  P^(k))  d|Jrk'  (U'2-10) 


for  almost  all  realizations  of  the  Markov  chain.   Thus  if  L(r^   ,  q    ) 

(k)  (k) 

is  averaged  over  r    produced  by  the  Markov  chain  and  q   '  generated 

by  the  straightforward  Monte  Carlo  scheme,  the  result  will  be  an 

approximation  for  G  . 

We  must  now  demonstrate  a  procedure  which  will  generate  a 
Markov  chain  which  satisfies  the  above  conditions.   It  will  be 
convenient  for  the  rest  of  this  section  to  deal  with  only  the 
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x -components,  x.,  of  each  3 "dimensional  vector,  r..   All  quantities 

(k) 
defined  above  in  terms  of  the  3k-dimensional  vector,  rv   = 

(r  ,  r  ,  . . .  ,  r  ),  r  =  r  =0,  are  assumed  to  have  analogous 
~1  ~2        ~k   ~0   ~k 

definitions  in  terms  of  the  k-dimensional  vector  x  =  (x,  ,  x„,  ...  ,  x  ), 

x„  =  x  =  Oo   The  generalization  to  three  space  components  will  be 

obvious,  and  the  notation  will  be  much  simplified  if  we  assume  all 

vectors  to  be  k-dimensional  and  drop  the  use  of  k  as  a  superscript 

or  subscript  to  indicate  dimensionality  for  the  rest  of  this  section. 

Let  w(x)  be  defined  by 


,k-l        ,k-l 


du     =  w(x)d       x,    d       x  =  dx,    dx_    . ..    dx.    , , 

x  ~  ~  ~  12  k-1 


(4,2.11) 


so  that 


w(x)   =   (2ir) 


k 
|2 


>2it 


exp 


S     Z      (x.-x.    .)' 
2  1      1-1 

1=1 


(4.2.12) 


A  transition  probability  function,    p(x' |x),    is   then  needed  to  satisfy 
Eq.    (4.2.8),   which  now  becomes 


/;•■/ 


P(x')  w(x')dk"1x' 


oo  oo 


.f ...   ff* .  . .  f  j>(x' \x)   P(x)w(x)  dk_1x  ■_  dk_1x'    . 


(4.2.13) 


-00  -00 
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First  define  the  function 


p*(x'|x)  =l±\ 


k-1 

2  exp 


-k 


k-1  /    x! 
Z  [x!  -  -1 

1=1  V1 


-1  +  xi+l 


(4.2.11+) 


This  function  satisfies  the  criteria  for  a  transition  probability- 
function  since 


p*(x' |x)  >  0 


(4.2.15) 


and 


00         00 


f'"f    P*(2'l;X)ak~V  =  !' 


-00        -00 


(4.2.16) 


The  function  p(x'  |x)'  is  now  defined  by 


P(x') 
p*(x«|x)   p^y-  ,  if  P(x')  <  P(x) 


P(x'|x)  -  <J  P*(x'|x)  +  5(x',x)  J  .  .  .      Jp*(x"|x) 

Cx"|P(x")  <  P(x)} 

if  P(x')  >  P(x)  , 


1  - 


P(x") 


pTxT 


,k-l  „ 
d   x  , 


(4.2.17) 


•here  Gx"|P(x")  <  P(x)}  is  the  set  of  all  x"  such  that  P(x")  <  P(x)  and 
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5  is  the  Dirac  delta  function.   The  function  p(x' |x)  is  a  valid 
transition  function  since  it  is  non-negative  -and  satisfies 

00  00 


I ...  rP(x,|x)dk"ix' 


-00  -00 


{x1 |P(x')   >  P(x)) 


p*(x' |x) 


+  B(x',x) 


[xM|P(x")  <  P( 


J^(-^)*k-v]^ 


r  p(x,)  ki 


{x'|P(x')  <  P(x)) 

=/   .  •  .       /W  te>*k~V 

(x'|P(x')   >  P(x)} 


f  ...    jrp*(x"ix)dk"ix*' , 

(x"|P(x")  <  P(x)) 


=  1. 


(4.2.18) 
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In  order  to  show  that  p(x' |x)  satisfies  Eq.  (4.2.13)  it  will 

T 
he  useful  to  define  the  "transpose",  x  ,  of  a  vector,  x  =  (x  ,  x  ,  ...,x  ), 

~  ~    •  l   d  k 

xQ  =  xk  =  0,  by  x±   =  xk_i,  i  =  0,  L,  .  ..,  k.   From  Eqs .  (4.1.29), 
(4.2.5),  and  (4.2.12)  it  can  be  seen  that 


w(x)  -  w(*  )  , 


(4.2.19) 


and 


P(x)  =  P(xX)   . 


(4.2.20) 


k-1       T    T  k-1  T 
Since  P(x)w(x)d   x  =  P(x  )w(x  )d   x  ,  and  since  the  integration  with 

respect  to  x  in  Eq.  (4.2.13)  is  carried  out  over  all  space,  it  is 

clear  that  Eq.  (4.2.13)  is  equivalent  to 


"       J    ' 


P(x,)w(x')dk"1x' 


■/•;•/ 


P(x')w(x*)    +   P(x'T)w(x,T)      ,     . 
~  ~  ~  ~  K.--L     , 

5 d         x.         > 


oo  oo 


■/■;•//-/ 


-00  -00 


p(x'  |x)   +  p(x'     |x    ) 


P(x)w(x)   dk_1x  dk_1x«       . 


(4.2.21) 
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As  a  lemma  for  the  proof  of  Eq.  (4.2.21)  it  is  shown  in  Appendix  A  that 


p*(x' |x)w(x)  =  P*(xT|x,T)w(x,T)  . 


(4.2.22) 


To  prove   that  p(x'|x)    satisfies  Eq.    (4.2.21)    it    is   sufficient  to   show 
that 


00  00 


v(x')   = 


-00  -00 


p(x' |x)   +  p(x'    |x    ) 


P(x)w(x)dk_1x    , 


=  2P(x')v(x')    , 


(4.2.23) 


where  v(x')  is  defined  by  the  above  equation.   Using  the  definition, 
Eq.  (4.2.17),  of  p(x' |x)  and  keeping  in  mind  that  P(x)  and  w(x)  are 
invariant  under  transposition  of  the  coordinates  of  x  it  follows  that 


v(x')  =f      .  .  .     J 
U|P(x)  >  P(x')} 


p*(x'|x)  +  p*(x'T|xT) 


P(x')w(x)dk'1x 


(x|P(x)  <  P(x')) 


p*(x' |x)  +  P*(x"T|xT) 


+  6(x',x) 


{x"|P(x")  <  P(x)} 


p*(x"|x)  +  p*(x"T|xT) 


1  -  pjy  I  dk""V  \    P(x)w(x)dk_1x  , 


(4.2.24) 
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where  the  integral  multiplying  the  delta  function  can  "be  extended  over 
the  set  {x"|P(x")  ^_P(x)}  since  the  integrand  vanishes  for  P(x")  =  P(x) 
Performing  the  integral  over  the  delta  function  v(x')  becomes 


v(x')  -  P(x') 


{x|P(x)  >  P(x')} 


p*(x'|x)  +  p*(x,T|xT) 


w(xjd   x 


{x|P(x)  <   P(x') 


p*(x*|x)  +  p*(x,T|xT) 


P(x)w(x)dk";Lx 


/ 


{x"|P(xM)  <   P(x')}   L 


p*(,x  x  )  +  p*(x   X   ) 


■/  •••  J 

{x"|P(x")  <   P(x')} 


P*(x"|x')  +  P*(x"T|x,T) 


P(x')v(x,)dk"1x' 


P(x")w(x')dk'1xn  . 
(4.2.25) 


Applying  Eq.  (4.2.22)  to  the  first  and  second  terms  in  the  last 
equation  one  finds  that  the  second  and  fourth  terms  cancel  and  there 


remains 


v(x')  =  P(x')w(x') 


(x|P(x)  >  P(x' )) 


p*(xT|x,T)  +  p*(xjx') 


d   x 


+  P(x*)w(x') J  ...       J 

{x"|P(x")  <  P(x'j} 


p*(x"|x')  +  P*(x'*T|xiT) 


,k-l  „ 

d   x 


(4.2.26) 


6o 


Since  the  sets  over  which  the  two  integrations  are  carried  out  are 
complements,  the  integrals  may  be  combined  into  an  integral  over  the 

whole  space.   Bearing  in  mind  that  the  variable  of  integration,  x,  may 

T 
just  as  well  be  replaced  by  x  one  obtains  the  result 


v(x')  =  P(x')w(x') 


00         00  00         00 

j ...  yV(x  ix'  )dk_ix  +  j . . .  yxp^(xTix,T)dk"ixT 


=  2P(x')w(x')  , 


(4.2.27) 


where  the  final  step  makes  use  of  the  normalization  condition  on 
P*(x|x'),  Eq.  (4.2.16), 

It  can  be  seen  from  the  definitions  of  p(x' |x)  and 
p*(x' |x)  that  there  can  be  only  one  ergodic  class  since  the  probability 

of  the  one  step  transition  from  any  vector,  x,  to  a  vector  in  the 

k-1 
volume  d   x'  about  x'  is  strictly  greater  than  zero  for  all  vectors, 

xJ,  satisfying  P(x')  >  0.   Equation  (4.2.10)  is  thus  satisfied  for  the 

p(x*  [x)  defined  above,  and  a  procedure  which  generates  a  Markov  chain 

with  this  transition  function  is  a  valid  importance  sampling  scheme 

for  the  present  problem.   The  transition  from  a  vector,  x,  of  the  chain 

to  the  next  vector,  x',  can  be  made  using  k-1  normally  distributed 

random  variables  and  one  uniformly  distributed  random  variable.   Let 


6l 


it    ii 
xo  =  xk 


o  , 


:i  =V^  Sj 


X  .  .  ,  +   X.  , 

1-1    1+1 


i  =  l,2,...,k-l  , 


(4.2.28) 


where  £.,  i=l,2, . . .,k-l,  are  normally  distributed  random  variables 
with  mean  0  and  variance  1.   If  P(x" )  <  P(x),  generate  a  random 
variable,  £,  which  is  uniformly  distributed  on  the  interval  (0,1)  and 
set 

P(x" ) 

(4.2.29) 


If  P(x")  >  P(x)  set 


~   '  "  *  ^  p(x)  ' 

x      ,    otherwise. 


(4.2.30) 


It  can  be  seen  that  the  probability  of  x"  given  x  is  P*(x" |x)  while 
the  probability  of  x1  given  x  is  p(x! |x). 
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4.3   Integration  over  the  Position  of  the  Third  Particle 

Besides  the  6(k-l)-fold  integration  arising  from  the 
approximation  to  the  Wiener  integral,  it  is  also  necessary  to  perform 
the  twofold  integration  over  the  position  of  the  third  particle. 
This  is  the  integration  over  (Pqj^q)  in  Eq.  (4.1.14),   Since  a  large 
number  of  Monte  Carlo  samples  are  necessary  for  reasonable  accuracy 
in  the  Wiener  integral  approximation,  it  is  impractical  to  compute 
G  (r  p,Q)  and  G  (r,p,Q)  at  each  point  of  a  2 -dimensional  grid  in  Q 
and  apply  a  standard  numerical  integration  technique.   If  we  perform 
this  last  twofold  integration  by  a  Monte  Carlo  technique,  however, 
the  sampling  may  be  combined  with  that  for  the  6(k-l)-fold  integrals. 
The  procedure  then  is  to  generate  M  samples  {(r;   ,  qj   ,  £.), 
t  =  1,2,...,M]  and  to  form 


M 

z 

t=l 


(ri2;M)  =  |  Z  L(r[k),  q[k),  Qt)  r"1^),  (l*.3.l) 


(k) 
where  L  is  defined  by  Eq.  (4.2.4),  r^   is  generated  from  the  Markov 

(k) 
chain,  q;.   is  produced  by  the  straightforward  sampling  scheme  for 

Wiener  paths,  and  £  =  (p_,0A,O)   is  picked  according  to  a  2-dimensional 

distribution 


r(g)d2S  -  rp(pQ)r0(eQ)dpQdeQ,  (M#2) 
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or 


r(S)a2£  =  rx(Qx)ry(Qy)d<yiQy,.  (4.3.3) 


in  Cartesian  coordinates.   If  Eq.  (4.3.3)  is  used  a  factor  of  - 
must  be  included  in  L  if  the  range  of  Q  is  -co  <  Q  <  <».   One  will 

y       y 

then  have  by  the  law  Of  large  numbers  that 


lim  Sd(r!2^M)  =  gd^r12^  (M-^) 

M-*oo 


with  probability  one,  where  g  (r ._)  represents  either  g  (r   )  or 
g  (r,p)  depending  on:  the  form  of  the  integrand  L(r   ,  q   ,Q  )• 

Including  the  integration  over  Q  in  the  Monte  Carlo  scheme 
will  increase  the  number  of  samples  necessary  for  a  given  accuracy 

since  L(r   ,  q   ,   Q)  will  have  an  increased  variation  on  being 

(k)      (k) 
taken  as  a  function  of  Q  as  well  as  of  r    and  q   .   The  variation 

with  respect  to  0^  turns  out  to  be  the  dominant  part  of  the  variance 

of  the  Monte  Carlo  scheme  at  moderate  temperatures  with  the  distri- 

(k)      (k) 
butions  of  r    and  £    being  fairly  sharply  peaked  about 

(k)    (k) 

r    =  q    =0.   This  variance  in  0  can  be  reduced  somewhat  by  making 

use  of  the  latitude  that  exists  in  the  choice  of  the  distribution 
r(§).   Qualitatively  the  variance  with  respect  to  Q  will  be  reduced 
if  the  global  behavior  of  T(Q)  is  similar  to  that  of  L(r/k  ,    q^    ',   £) 


6k 


as  a  function  of  Q.   Since  L  is  also  a  function  of  r    and  q 
computing  time  restrictions  indicate  that  it  would  probably  not  be 
worthwhile  to  develop  an  elaborate  importance  sampling  scheme  in  the 
variable  Q.   The  choice  of  r(Q)  will  therefore  be  restricted  to 
distributions  which  can  be  generated  easily  from  uniform  or  normally 

distributed  pseudo  random  numbers. 

(k)   (k) 
In  order  to  get  an  idea  of  the  form  of  L(r   ,   q   ,  Q) 

function  of  £  one  must  consider  g  (r,p)  and  g  (r,  p)  separately.   For 

the  function  g  (r,p)  the  integrand  is 

_  t  (k)  (k)  ns   r  j^3  (k)  x  (k)\   n 

Lc(£    i  a    '  s)  =  pi3 r-fi"    +  |r      -  1\ 


as   a 


^3(k)-h(k))-i- 


(i».3.5) 


while  for  g   (r,p)   the   integrand   is 


V^S(^S)=[Dl3(^a(k)^£(k))-Di3 


(,(k)) 


D       KiQW        I  rU)\    -    D        (q(k))l 
23  I    2    S  2  ~        j  23    VS        ; 


(^.3.6) 
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*.v  *<    4-  -v  +•     *      (k)   j   (k)        i  i       (k)    (k) 
Since  the  distributions  of  r    and  q    are  peaked  near  r    =  q    =  0, 

the  function 


Lc(0,0,Q)  = 


exp  {-PV(Q  +  ^)  -  |^  V"(Q  +  ^f-)}    -  1 


X 


exp  {-PV(Q 


'12. 


gx: 


-12 


2nk   ^    2  n 


(4.3-7) 


gives  some  idea  of  the  variation  of  L  with  respect  to  Q.  Equation 
(4.3*7)  is  just  the  second  order  Wigner-Kirkwood  approximation  to  the 
classical  cluster  function 


f13f23  =  (e    13J  -  l)(e    23   -  1)   , 


(4.3.8) 


so  that  the  behavior  of  this  classical  function  with  respect  to  Q  will 

give  some  idea  of  the  global  structure  of  L  as  a  function  of  Q. 

c  ~ 

Figure  2  shows  a  rough  contour  map  of  f  ^fpo  for  T  =  20  K 
and  r  p  =  1.1a,  where  a  is  defined  by  Eq.  (2.1.2).   The  structure  of 
this  cluster  function  is  too  complex  to  attempt  to  fit  it  in  detail 
with  a  simple  probability  distribution.   It  can  be  seen,  however, 
that  the  structure  of  the  function  is  centered  about  Q  =  0  and  that 
the  function  becomes  negligible  at  a  distance  of  2  or  3  times  a  in 
the  Q  direction  and  at  a  slightly  larger  distance  in  the  Q 


direction."'  A  reasonable  sampling  distribution  for  such  a  function 
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FIGURE  2.  fi3(ri2iQ)Xf23(r12i9)  AS  A  FUNCTION 
OF  Q  FOR  r12  =  l.lcr "AND  T=20°K. 
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might  be  a  2 -dimensional  normal  distribution  with  a  standard  deviation 

on  the  order  of  a   in  the  Q  direction  and  a  standard  deviation  in  the 

T 

Q  direction  on  the  order  of  a   +  const,  x  rnr,.   The  distribution  used 
x  12 

for  the  integration  over  Q  in  the  numerical  computations  to  follow  was 


T^)i%  "^7xexp[  k)  A  exP  I"  *  J dVQy  '     (M'9) 


where  the  standard  deviations  7  and  7  were  chosen  empirically  to 
minimize  the  variance  of  g  (r  ;M)  for  a  small  value  of  M. 

If  one  tries  to  take  r'1^  =  q^k'  h  0  in  L  (r^\    q_^\    Q) 

~      ^,  e  ~     ~     ~ 

-3 
the  function  vanishes.   In  fact  this  function  vanishes  to  order  T 

in  the  Wigner-Kirkwood  expansion,  Eq.  (3.2.25)^  and  it  is  thus  to  be 

expected  that  the  correction  term  g  (r,p)  will  be  negligible  except 

at  low  temperatures  and  for  small  r,p.   In  order  to  find  an  efficient 

sampling  scheme  for  g  (*\p)  advantage  can  be  taken  of  several  symmetries 

of  the  integrand.   The  function  g  can  be  written  as  the  sum  of  two 

integrals 


ge(r12)  .  £\Il2)   ♦  g(2'(r12)   ,  (U.3.10) 
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where 


41)<ri2)  ■  crfer  M 

d.    ~  ~ 

x  E  JF12(r(T))  |F13(^  4(t)  +  ^r]  -  F13(4(t))|  |r(l)  =  3(1)  =  o}, 


['«(' 


iF12(r(T))  I  F23  f  ^  q(-r)  -  :L-^-\  -   F^qCx)) 


r(l)  =  q(l)  =  Oj, 

(^3.11) 


and  where  F  ,  F,_,  and  F   are  defined  by  Eq.  (4.1.11).   Both  g'  ' 

(2) 

and  g    are  functions  of  the  magnitude  r  p  only,  so  r,p  may  be  replaced 

(2) 
by  -  r,p  wherever  it  appears  in  g   (r,p).   Furthermore,  Wiener  measure 

on  r(x)  is  invariant  to  the  replacement  of  r(r)  by  -r(x).   If  one  makes 
the  two  changes  of  variables  r,p-»  ~£-io  an(^  r(T)~>  ~£(T)  ^n  §   >  then 
F,p(r(T))  is  unchanged  due  to  the  spherical  symmetry  of  V(r)  while 

F23  ~2  S(t)  "  Z~2~y   Fl3  (  ~2  S(t)  +  Z~2~ 


and 


F23(q(T))-Fi3(c/T)) 
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Thus 


^(V  =  g'2^) 


(It.  3.12) 


and  g  (rip)  may  be  written  as 


(r10)  - 


■  12'   W2(l,2) 


d3Q 


X  E^F12(r(T)) 


(t) 


Fl3l~2^T)  +: 


-  F13(q(T)) 


r(l)  =  q(l)  =  O; 

(^•3.1 


In  this  last  form  a  more  natural  origin  for  the  integration 


*12 


over  Q  is'Q=-  — ==-,  so  let  the  variable  of  integration  he  changed  to 


Q'"  =  Q  + 


-12 


(4.3.1*0 


Then 


ge(r!2)  =  J  dV  Ge(r12'S!)   * 


(4.3.15) 


where 


TO 


Ge(r12,£')  -  2E 


1i 


F22^r^ 


W2(l,2) 


F13[  Q',  -|c 


p(t) 


Fi3(S'^(t)) 


and  where 


r(l)  =  q(l)  .  (A  , 


(^•3.16) 


F13(s;a(^))  -  exp 


P  /  VI7-  q(T)  +  Q' 

.oJ   V*  ~ 


(t.3.17) 


The  prime  on  Q'  will  be  dropped  for  the  remainder  of  the  discussion 

of  g  (r   ).   In  spherical  coordinates  the  cylindrical  symmetry  in 

cpn  is  not  altered  by  the  change  in  origin  provided  the  polar  axis 

remains  in  the  Q  direction, 
x 

The  integrand  of  g  (r  _)  is  the  difference  of  two  terms 
depending  on  g.   This  difference  will  be  small  if  both  terms  are 
small  or  nearly  equal  so  we  consider  the  behavior  of 


F13(Q,0)  =  e 


■PV(Q) 


(4.3.18) 


with  respect  to  Q.   This  function  is  the  limit  of  both  of  the  terms  in 


the  square  brackets  in  Eq.  (4.3-16)  for  r(f)  =  q(i)  =  0.   For  pn  <  a, 
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F, _(Q,0)  goes  to  zero  rapidly.   The  function  has  a  maximum  at  p  -  1.12a 

'    6  \  ^ 

and  approaches  one  asymptotically  as  I  1  +  — z"  1  for  p  »  a.  The 

*V 

integrand  in  Eq.  (4.3-15)  should  thus  be  zero  for  pn  «  o  and  approach 
zero  again  for  pn  »  a  as  both  terms  in  the  difference  approach  unity. 
A  simple  probability  distribution  with  this  same  qualitative  behavior 

is  the  x  distribution  with  n  degrees  of  freedom  for  n  >  k.      For  the 

2 
£  integration  in  g  (rip)  then  a  x  distribution  with  12  degrees  of 

freedom  was  used  for  pn  and  £~   =   cos  9n   was  chosen  from  a  uniform 

distribution,  -1  <  £n  <   1.  Thus  we  define 

re(S)d2Q  =  ^2  exp  I"  -H  dpQd|Q'         (^3*19) 


The  number  of  degrees-  of  freedom  and  the  value  of  y     were  chosen 

P 

empirically  during  preliminary  computations  to  minimize  the  variance 
of  g  (r  p;M)  for  small  M. 

In  initial  calculations  of  g  (r1p;M)  it  was  found  that 
g  (r  p)  was  indeed  small  compared  to  g  (r   )  but  that  the  variance  of 
g  (r  p;M)  was  quite  large.   It  was  thus  necessary  to  apply  a  further 
variance  reducing  technique  to  the  Monte  Carlo  sampling  in  the 
calculation  of  g  .   The  technique  used  was  the  method  of  antithetic 
variates  described  by  Hammersly  and  Handscomb  [26],  and  although  it 
applies  to  the  Monte  Carlo  sampling  as  a  whole,  it  is  best  described  in 
connection  with  the  sampling  in  the  Q  integration.   The  principle  of 
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the  antithetic  variates  method  is  to  generate  several  correlated 
samples  simultaneously  in  such  a  way  that  their  sum  has  a  considerably 
smaller  variance  than  the  samples  taken  independently.   In  the  present 
case  it  is  recognized  that  g  (r  „)  is  small  so  that  one  tries  to, 
correlate  the  samples  in  such  a  way  that  their  sum  is  near  zero. 
From  Eq.  (4.3-16)  it  can  be  seen  that  the  cylindrical 
symmetry  of  G  (r 1P>Q)  about  the  Q  axis  requires  that 

Ge(r12'-S)  =  Ge(r12'pQ'*Q  +f'V'  {k'3-20) 


Thus  the  range  of  6Q   can  be  restricted  to  0  <  9Q  <  .ny   (0  <  £Q  <  l) 

and  the  integrand  taken  as  the  sum  of  G  (r,p,Q)  and  G  (r  „,-<^). 
Another  pairing  of  samples  can  be  made  due  to  the  fact  that  the 

probability  of  a  Wiener  path  q(x)  is  the  same  as  the  probability  of 

(k) 
-q(-T).   Since  successive  path  approximations  q    are  chosen 

(k)       (k) 
independently  the  samples  for  q    and  -q    may  be  combined. 

From  Eqs,  (4.1.28),  (4,1.29),  and  (4.2,4)  it  is  seen  that 

the  integrand  of  the  multiple  integral  approximation  for  G  (r.p,£) 

as  given  in  Eq.  (4.3.16)  is 


Le(£w,  s«  g)  =  2  ^13(S.  f  a(k)  +  -2-)-  \^>  s(k)]  > 
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where 


D13(^S(k)) 


=  exp 


*i 


1=1 


{r«^)+kY[r« 


./^+5i-i 


h-Q 


^ltH  +  =Yl/. 


x  3i+Si-n 


+  Q 


(4.3.22) 


(k)  (k) 

Then  if  r    is  generated  from  the  Markov  chain,  q)        is  generated  by 

the  straightforward  scheme  for  Wiener  paths,  and  Q,  is  chosen  from  the 

distribution  r  (Q)  defined  by  Eq.  (4.3.19)  but  with  |n  restricted  to 

the  range  0  <  £_  <  1,  then  with  probability  one 

ge(r12)  -llB-JZ  A(r[k),  q[k\  Qt)  r"1^)  ,  (4.3.23) 

M-*»    t=l 


where  the  antithetic  variates  sample 


./  (k)   (k)  ^N 
A(£   >    <1   ;  S) 


is  defined  by 


V£(k).  S(k)-  S)  -  W°.  a'"''  -S) 


A(r<*\  *«  S)  =  | 


(U. 3.2"0 


7^ 


This  sum  should  be  more  nearly  zero  than  its  component  terms  since 
the  portion  of  L  which  is  antisymmetric  in  Q  will  cancel  between 
terms  1  and  2  and  between  terms  3  and  k   on  the  right  side  of  Eq. 

(k.3.2k)   while  similar  cancellation  will  take  place  for  the  component 

(k) 

of  L  which  is  antisymmetric  in  the  variable  q    between  terms 

1  and  3  and  between  terms  2  and  h. 
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5.    NUMERICAL  RESULTS 

5 .1   Results  and  Accuracy  Considerations 

Using  the  methods  developed  in  the  preceding  sections 
g  (*%„;  20,000)  and  g  (r  ',    20,000)  were  calculated  for  temperatures 
T  =  5°K,  10°K,  20°K,  30°K,  40°K,  50°K,  75°K,  100°K,  and  273.l8°K 
and  for  values  of  r   in  the  range  0.6  <  r.p/a  <  4.0.   The  third 
virial  coefficient  C(T)  was  computed  for  all  "but  the  lowest  temperature, 
5  K,  using  Eq.  (2.3°20)°   A  straightforward  numerical  integration 
was  performed  over  r,p  with  g, (r^p)  approximated  by  g  (r,p;  20,000)  + 
g  (r,p;  20,000)  and  with  the  values  of  Wp(r,p)  and  bp  given  in  [6], 
The  computations  were  performed  on  the  ILLIAC  II  computer  located  at 
the  University  of  Illinois.   The  details  of  the  computation  and  a 
comprehensive  table  of  results  are  presented  in  Appendix  Be   A  summary 
of  the  results  for  g-^r^;  20,000)  =  gc(r12;  20,000)  +  g  (r12;  20,000) 
over  the  range  of  temperature  is  shown  in  Figure  3>  and  the  computed 
third  virial  coefficient  C(T)  is  plotted  versus  temperature  in  Figure  h< 

i 

There  are  two  sources  of  computational  error  which  cause 
g  (rip;  20,000)  and  C(T)  to  differ  from  the  correct  theoretical  values 
of  g-,(r,p)  and  the  third  virial  coefficient  for  the  interaction 
potential  of  Eqs .  (2.1.1)  and  (2,1,2).   These  are  the  error  introduced 

by  the  k-fold  integral  approximation  to  the  Wiener  integrals  and  the 

(k) 

Monte  Carlo  sampling  error.   The  order  k  of  the  approximations  r 

(k) 

and  q    to  the  Wiener  paths  r(t)  and  q(f)  was  chosen  for  each 

temperature  T  so  that  Tk  >  260  K,   Thus  the  accuracy  of  the  multiple 
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FIGURE  3.     X    g^hz  ; 20,000)  AS  A  FUNCTION  OF  r^/C  AND  T, 
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integral  approximations  at  any  temperature  is  at  least  equal  to  the 
accuracy  of  the  two  term  Wigner-Kirkwood  approximation  (first  two  terms 
in  Eq.  (3 • 3<25))  at  260  K.   As  an  estimate  of  this  accuracy  the  third 


term  in  the  Wigner-Kirkwood  approximation  to  the  second  virial 
coefficient  of  He  at  26o°K  amounts  to  only  0.7$  of  the  total.   It 
will  be  seen  that  an  error  of  this  magnitude  is  negligible  compared 
to  the  Monte  Carlo  sampling  error. 

As  an  estimate  of  the  Monte  Carlo  sampling  error  the  sample 
standard  deviations  of  the  sample  means,  g  (r,p;  20,000)  and 
g  (r,p;  20,000),  were  computed.   The  numerical  results  for  g,  (r-.p;  20,000) 
and  the  associated  error  estimates  are  shown  plotted  against  r,p/a 
for  T  =  5°K,  20°K,  75°K,  and  273<l8°K  in  Figures  5,  6,  7,  and  8, 
respectively.   The  complete  set  of  graphs  can  be  found  in  Appendix  B. 
The  length  of  the  error  bar  associated  with  A  g  (r,p;  20,000)  is 
equal  to  twice  the  sum  of  the  sample  standard  deviations  of 
^"3  g  (r,0j  20,000)  and  \"3  g  (r10j  20,000), 

C   -Li:  e   ±C. 

Although  the  Monte  Carlo  sampling  connected  with  the  Wiener 
integration  and  with  the  Riemann  integration  of  G(r  p,  Q,)  over  the 
position  Q  of  the  third  particle  is  done  as  a  single  process,  it  is 
convenient  to  consider  the  sampling  error  to  be  a  sum  of  separate 
contributions  from  these  two  sources.   The  sampling  error  in  the 
integration  with  respect  to  Q  depends  on  the  variation  of 
G(r,p,  Q)  T   (Q)  as  a  function  of  Q.   Since  G(r,p,  Q)  contains  functions 
of  the  form  e    ^~^12'  '   which  have  an  increasing  variation  in  Q,  as  |3 
increases,  it  is  to  be  expected  that  the  sampling  error  due  to  the  £ 
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integration  will  increase  with  decreasing  temperature. 

A  factor  of  (3  also  multiplies  the  exponent  of  the  functional 
integrand  of  the  Wiener  integration  so  that  the  variance  of  the  Monte 
Carlo  procedure  for  computing  the  multiple  integral  approximation  to 
the  Wiener  integral  will  also  increase  as  T  decreases  as  a  result  of 
the  increase  in  (3.   The  variance  of  the  multiple  integral  approximation 
is  further  influenced  by  the  fact  that  the  Wiener  random  functions 
r(  t)  and  q(l)  in  Eqs .  (4.1,8)  and  (4ol.ll)  are  multiplied  by  the  thermal 
wavelength  \.      Since  \  increases  with  decreasing  temperature  the  variance 
again  increases  as  T  decreases.  A  final  contribution  to  the  same 
behavior  of  increasing  variance  with  decreasing  T  results  from  the  fact 
noted  above  that  k  was  increased  with  decreasing  T  in  order  to  maintain 
the  accuracy  of  the  multiple  integral  approximation  to  the  Wiener 
integrals.   The  multiplicity,  6(k-l),  of  the  integral  thus  increases  with 
decreasing  T  resulting  in  an  increasing  variance  for  the  Monte  Carlo 
procedure. 

The  standard  deviation  of  the  multiple  integral  results  also 

(k) 

depends  upon  the  choice  of  the  initial  vector  r     of  the  Markov 

chain  {r,  '   )«   Due  to  the  one  step  memory  inherent  in  the  chain  the 

(k) 
choice  of  an  initial  vector  r  '   which  is  improbable  with  respect  to 

the  probability  density  P(r^k')du  v  (Eqs.  (h.1.26)   and  (U.2.5))can 

(k)  ' 
cause  a  large  variation  in  the  integrand  as  the  vectors  r     move 

(k)  (k) 

from  r     toward  the  region  of  maximum  P(r'   )du  ,  .  An  attempt  was 
~o  °  ~   '      rk 

(k) 
made  to  reduce  the  error  caused  by  a  bad  choice  of  r     for  each  pair 

J  ~o  * 

(k) 

of  values  of  r,Q  and  T.  At  high  temperatures  P(r   )du  ,  assigns  the 


81+ 


(k) 
highest  probability  to  vectors  in  the  neighborhood  of  r    =  0  provided 

that  P(r    )  is  slowly  varying  over  the  range  of  r,p  +  r    .   Thus  for 


(k) 
be  r     =  0.   At  lower  temperatures  the  initial  vector  for  the  largest 
~o  ° 


the  highest  temperature  and  largest  r,p  the  initial  vector  was  taken  t< 

l1  vi 

(k) 

value  of  r,  ~  was  taken  to  be  the  last  vector  r.„    generated  by  the  Markov 
12  ~M    D         J 

chain  for  the  same  value  of  r   at  the  next  higher  temperature.   For 
successively  smaller  values  of  r,p  the  initial  vector  was  taken  to  be 
the  last  vector  generated  at  the  same  temperature  and  the  next  higher 
value  of  r   . 

It  can  be  seen  from  the  tables  in  Appendix  B  that  g  (r,p;  20,000) 
is  small  with  respect  to  g  (r,p;  20,000)  as  was  to  be  expected  from 
the  fact  that  g  (r,  p)  vanishes  to  order  T   in  the  Wigner-Kirkwood 
expansion,  \     g  (r, p;  20,000)  is,  on  the  average,  an  order  of  magnitude 
smaller  than  X     g  (r,p;  20,000).   The  standard  deviations  of  these 
quantities,  however,  are  of  approximately  the  same  size  so  that  the 
effect  of  adding  g  to  g  in  the  approximation  for  g,  is  primarily  an 
increase  in  the  expected  error  in  g,  . 

The  quantities  g  (r, p;  M)  and  g  (r  p;  M)  were  also  computed 
with  M  <  20,000  at  several  pairs  of  values  for  r,p  and  T,   The  standard 
deviations  in  these  quantities  were  observed  to  decrease  as  M  2  as 
would  be  expected  of  a  well -formulated  Monte  Carlo  procedure. 
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5.2   Numerical  Comparisons  with  Theory  and  Experiment 

It  can  be  seen  from  Figure  3  that  the  position  of  the 
minimum  in  g ,(r,p)  remains  fairly  constant  at  r   -  l„8a  over  the 
entire  temperature  range .   The  position  of  the  minimum  in  g  (t,0) 
can  "be  interpreted  in  terms  of  a  crude  argument  about  the  spatial 
arrangement  of  three  particles .   If  one  takes  the  most  probable 
arrangement  of  particles  as  a  linear  arrangement  with  each  particle 
falling  at  the  minimum  of  the  potentials  due  to  its  two  nearest 
neighbors,  then  for  the  Lennard-Jones  (12^6)  potential  there  should 
be  maximum  probability  of  finding  a  particle  at  distances  of  1.12a 
and  2.2ka   from  a  reference  particle.   Correspondingly  there  should 
be  a  minimum  probability  of  finding  a  particle  midway  between  these 
two  positions  or  at  1.69a.   If  Wp(r,p)  is  slowly  varying  near  this 
minimum  of  the  pair  distribution  function,  then  the  corresponding 
minimum  in  g  (r, p)  should  fall  at  approximately  the  same  value  of 
r,p.   The  calculated  position  of  the  minimum,  1.8a,  agrees  well 
enough  with  this  result  to  indicate  that  the  minimum  can  be 
interpreted  physically  as  the  gap  between  the  first  and  second  shells 
of  atoms  about  the  reference  atom. 

Classical  and  two  term  Wigner-Kirkwood  approximations  for 
g  (r   )  have  also  been  claculated  using  respectively  the  first 

-pa, 

and  first  two  terms  of  Eq.  (3°3«25)  for  <1,2,3  |e   J |  1,2, 3>. 


86 

The  integration  over  the  position  of  the  third  particle  was  carried 
out  by  a  straightforward  iterated  Simpson's  rule  approximation. 
Comparisons  between  the  classical,  Wigner-Kirkwood,  and  quantum 
calculations  are  shown  in  Figures  9,    10,  11,  and  12.   The  classical 
values  of  g,  (r  p)  have  previously  been  calculated  by  Henderson  [28], 
and  the  results  presented  agree  with  his  calculations  wherever  there 
is  overlap.   In  observing  these  comparisons  it  should  be  noted  that 
when  g  (r   )  is  used  to  calculate  a  physically  measurable  quantity 
it  is  multiplied  by  Wp(r,p)  which  goes  to  zero  rapidly  with  decreasing 
r  p  for  r.p/a  <  !•   Thus  values  of  g  (r.p)  for  r,p  <  0.8a  are  not 
significant. 

The  two  term  Wigner-Kirkwood  approximation  is  fairly 
accurate  down  to  about  T  =  10  K.  The  position  of  the  minimum  in 
g,(r.p)  corresponding  to  the  intershell  gap  remains  approximately 
constant  for  the  Wigner-Kirkwood  data  while  a  significant  shift 
toward  smaller  values  of  r,p  with  decreasing  temperature  occurs  in 
the  classical  position  of  this  minimum.  At  5  K  the  Wiener  integral 
method  gives  results  which  differ  markedly  from  the  two  term 
Wigner-Kirkwood  approximation.   Results  were  obtained  for  only  a  few 
values  of  r,p  at  T  =  5  K  due  to  the  amount  of  computation  time 
involved  in  the  Monte  Carlo  procedure  for  large  values  of  the  order, 
k,  of  the  Wiener  integral  approximation  (k  =  52  at  T  =  5  K; . 
Considering  the  a  posteriori  information  on  the  accuracy  of  the 
two  term  Wigner-Kirkwood  approximation  it  would  seem  that  reasonably 
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accurate  results  could  be  obtained  with  Tk  >  30  K  rather  than 

Tk  >  260  K  so  that  it  would  be  possible  to  repeat  the  5  K  computation 

in  one  tenth  the  time. 

There  seems  to  be  no  experimental  data  on  the  pair 

k 
distribution  function  in  He   gas  at  low  temperatures.   The  pair 

distribution  function  for  liquid  helium,  however,  has  been  measured 

by  neutron  diffraction  (Hens haw  [29])  for  several  values  of  the 

density.   Figures  13  and  lU  show 


W2(r12)  [1  +  n  gl(r12)]  -  d  ^     =  — &-  (5-2.1) 

n  o 


for  T  =  5  K  and  number  densities  n  corresponding  to  mass  densities 
p  =  0.095  gm/cc  and  p  =  0.18*4-  gm/cc  respectively  compared  with 
Henshaw's  data  for  these  same  densities.   Since  only  3 -particle 
effects  have  been  taken  into  account  in  the  computation,  agreement 
cannot  be  expected  beyond  the  second  peak  in  the  distribution  function, 
To  indicate  the  variation  of  p(r,„)/p  with  temperature,  smooth 
curves  have  been  fitted  to  the  numerical  results  for 


W2(ri2)  [1  +  n  g^r^)]  -  -j±-  (5-2.2) 


for  temperatures  10  K,  20  K,  and  i+0°K  at  a  density  of  pQ  =  0.184  gm/cc. 
The  resulting  curves  are  shown  in  Figure  15 • 
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Two  factors  influence  the  comparison  between  Henshaw's 
neutron  diffraction  data  for  liquid  He  and  the  numerical  results 
of  this  work.   The  first  is  the  choice  of  the  Lennard-Jones  (12-6) 
potential,  Eq.  (2.1.2),  as  the  correct  intermolecular  potential  for 
helium.   The  second  is  the  validity  of  the  cluster  expansion  for 
temperatures  and  densities  at  which  helium  is  a  liquid*   The  influence 
due  to  each  of  these  factors  cannot  be  separated  without  recomputing 
Wp(r  0)   and  g  (r,p)  for  a  different  intermolecular  potential.   However, 
barring  compensating  errors  in  these  two  approximations,  Figures  13 

and  ~\.K   indicate  that  both  approximations  are  fairly  good  for  the  pair 

k  o 

distribution  of  liquid  He  at  5  K. 

Using  Eq.  (2.3-20)  the  third  virial  coefficient  C(T)  was 
evaluated  for  all  temperatures  at  which  g-,(r,p)  was  calculated  except 
for  5  K.   C(5  K)  was  not  calculated  partly  due  to  the  large  sampling 
errors  in  g  (r   )  at  5  K  but  primarily  due  to  the  fact  that  values 
of  g  (r   )  for  r   >  2 =8 a  give  a  large  contribution  to  the  integral 
for  C(5  K)  and  no  numerical  results  were  obtained  for  this  range 
of  r1C)«   Figure  l6  shows  the  numerical  results  for  C(T)  compared  to 
the  values  adopted  by  Keesom  [30]  for  the  experimental  third  virial 
coefficient  and  the  values  calculated  classically  from  the  Lennard-Jones 
(12-6)  potential  (Hirschf elder  [1]). 

Over  almost  all  of  the  temperature  range  the  computed  values 
fall  consistently  below  the  experimental  results.   Since  Eq.  (2.3>20) 
does  not  depend  on  the  neglect  of  higher  order  terms  in  the  cluster 


96 


rf_cmL\ 
4  uVmole2/ 


600- 


500- 


400  ■■ 


300  ■ 


200  - 


100  ■ 


EXPERIMENTAL  DATA  ,  KEESOM    [30] 


THIS   WORK 


CLASSICAL    RESULTS  ,  HIRSCHFELDER   [l] 


50 


+— *T(°K) 


100  150  200  250  300 


FIGURE   16.   COMPARISON  OF  CALCULATED  THIRD  VIRIAL  COEFFICIENT 
WITH   CLASSICAL   AND  EXPERIMENTAL  VALUES. 


97 

expansion  the  cause  of  this  difference  is  evidently  the  choice  of  the 
pairwise  additive  Lennard- Jones  (12-6)  potential  as  the  correct 
3-particle  potential  for  helium.   Some  recent  calculations  by  Sherwood, 
De  Rocco,  and  Mason  [31]  indicate  that  nonadditive  3-body  forces 
have  a  sizeable  effect  on  the  third  virial  coefficient.   The  ways  in 
which  g  (r   )  and  C(T)  differ  from  the  experimental  values  can  be 
reconciled  qualitatively.   Figures  13  and  Ik   indicate  that  the 
calculated  g, (r,p)  is  fairly  consistently  larger  than  the  experimental 
g  (r  p).   On  performing  the  integration  over  r   in  Eq«  (2.3»20) 
this  would  lead  to  a  C(T)  which  would  be  smaller  than  its  experimental 
values  as  is  verified  in  Figure  l6„ 

5«3   Conclusion 

The  comparisons  of  Section  5=2  point  to  two  conclusions 

k 
regarding  the  physics  of  the  He   system.   The  first  is  that  quantum 

mechanical  calculations  based  on  the  cluster  expansion  and  assuming 

a  pairwise  additive  Lennard -Jones  (12-6)  potential  can  produce  a 

reasonable  qualitative  fit  to  the  pair  distribution  function  in 

k 
liquid  He  .   Secondly,  the  consistent  deviation  of  the  calculated 

third  virial  coefficient  from  the  experimental  values  over  the  full 

temperature  range  indicates  that  the  choice  of  the  latter  potential 

k 
for  the  3-body  interaction  in  He  is  not  good  enough  to  yield 

quantitative  agreement  with  experiment, 
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From  the  mathematical  viewpoint  it  has  "been  shown  that  the 
Wiener  integral  may  "be  a  useful  computational  tool.   The  upper  bounds 
of  Section  3-2  indicate  that  the  expression  of  a  quantity  in  terms 
of  Wiener  integrals  can  yield  qualitative  and  semiquantitative 
information  on  the  size  and  behavior  of  that  quantity.   Finally,  the 
numerical  evaluation  of  a  Wiener  integral  expression  has  been  shown 
to  be  computationally  feasible  in  a  situation  where  the  amount  of 
computation  involved  in  integrating  the  original  partial  differential 
equation  giving  rise  to  the  Wiener  integral  would  be  prohibitive. 

With  the  numerical  techniques  developed  it  would  be 
interesting  to  investigate  the  effect  on  the  third  virial  coefficient 
of  a  change  in  the  intermolecular  potential.   A  quantum  mechanical 
calculation  of  the  third  virial  coefficient  might  offer  a  way  of 
comparing  the  potentials  which  are  indistinguishable  with  respect  to 
the  second  virial  coefficient.  Another  possible  extension  of  this 
work  would  be  an  investigation  of  the  behavior  of  the  exchange  terms 

in  the  3-particle  density  matrix.  Of  course,  both  direct  and 
exchange  terms  can  be  computed  for  systems  of  more  than  three  particles 
using  the  same  formulation,  but  indications  are  that  the  necessary 
computation  increases  greatly  with  each  added  particle. 
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APPENDIX  A 


LEMMA  FOR  THE  PROOF  OF  Eq.  (4.2.21) 
Lemma:   p  (x'|x)w(x)  =  p  (x  |x'  )v(x'  ) 


(4.2.20) 


Proof:   From  Eqs .  (4.2.11)  and  (4.2.13)  it  follows  that 


p  (x  x   jw^x  ) 


k-1 


k 

2 


exp 


k-1 
-k   Z 
1=1 


T       ,T 
x.  .  +  x.  . 

i-l    i+l 


i=l  J 


(A.l) 


T   T 
By  the  definition  of  x  ,  x.  =  x   .,  i=0,l, ...,k,  so  Eq.  (A.l)  becomes 

~     1     K.- 1 


p  (x  x   )w(x'  ) 


i 


k 
2 


exp 


[-k  1  (vi 


x.  .  ,  +  xn  .  ,\p 

k-i+1    k-i-1*^ 


i=l  -1 


(A.2) 
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Setting  £   =  k-i,  one  finds  that 


p  (x  |x'  ;w(x'  ; 


k\   2 


k-1     1 
2 


(2«r 


27'2 


k-1 


k-1 


-k   L  Xo  +  k   L  x  pX  „   -, 
£=1     l  £=1       l  l 


k-1 
k   Z 
£=2 


x  „x 


r£-± 


k  k£x  ili+i +  xi-i} 


? 


i=i 


h 


^  i=o   * 


+  k  Jli     X£X£+1   "  2  f=1      ^ 


(A. 3) 


and  making  use  of  the  fact  that  x  =  x,  =  x1  =  x'  =  0,  this  "becomes 
&  o    k    o    k 


p  (x  x'  Jw^x  ) 


k-1     1 

k^2  '  to    \2  l*\2 

(2k)  [^)    exp 


k-1  k-1 

-k   Z   (xjr  +  k   Z   x'. 


L    i=l 


j=l 


J  j-l 


k-1 
+  k   Z   x'.x .  , 

.1=1   J  J+1 


k  kz'    (x^i +  xi-iV 


i=l 


k-1   2     k 
-k   Z   x  „  +       Zx.x. 


i=l        J=l 


J  j-l 


(AA) 


where  j  =  i+1.   Equating  the  dummy  subscripts  and  combining  the  sums, 
one  obtains 


10^ 


P  (x  |x'  )w(x'  ) 

k-1     1   /  \  k 


^•2  (-)2^)2exp[-k|fx.-!!M4ikX2 


2   ^   (xi  "  Xi-1)  J   ' 


p  (x1 |x)w(x)  ,  (A. 5) 


which  completes  the  proof, 
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..APPENDIX  B 


DETAILED  COMPUTATIONAL  PROCEDURE  AND  NUMERICAL  RESULTS 
FOR  gc(r12;M),  ge(r12;M)  and  C(T) 


The  step  by  step  procedure  for  computing  the  approximation 
g^(r12;M)  to  g^(r^)  is  as  follows: 

1.   Generate  M  independent  Wiener  path  approximations 

(k) 
consisting  of  3k-dimensional  vectors  q)_    ' 3    t  =  1,2,.,.,M,  drawn 

from  a  probability  distribution  du   by  Eq.  (^.2.2),  q  =  q  =  0, 


3i 


'li1  -  0 


1  - 


i-i  +ii 

k 


W\ 


i  -  1 


,  i  =  L,2,...,k-1   , 


(B.l) 


where  the  |.  are  3_dimensional  random  variables,  the  coordinate 
variables  of  which  are  independent  and  normally  distributed  with 

mean  0  and  variance  1. 

(k) 
2.   Generate  an  M  step  Markov  chain  [r)_        t=l,2,  .  . .  ,M) 

having  the  stationary  probability 


D12(£(k)) 
WQ(1,2)    d^rk 


by  first  generating  r;    using  the  procedure  described  by  Eq,  (^.2.28) 


io6 


ind  setting  r^*'  =  r^    if  D^r^')  >  D^fet-l^  while  if 


D12^tk))<D12^t-|)Set 


(k)        "(k)      •+>,         v,k-t+      °12(£t      ) 
r>        =  iv        with  probability  7TY~ 

D      (r{    }) 


and 


D     (r(k)) 
(k)           (k)      ...           KM1„      ,        °12^t      } 
rl        =  r :        with  probability   1   -  ,    \ 

D   (rl k)) 
12v~t-ly 


3-   Generate  M  independent  random  vectors  Q  =  (Q  ,Q  , 0)  , 

~t     x  y   t 

2 

t  =  1,2, ... ,M, according  to  the  probability  T   (Q)d  Q. 


h.      Form  the  average 


M 
(^JM)  ..  i   £  Lc(r«  ,  q<k»  ,  Qt)  r;1^)  .         (B.2) 


The  equations  pertinent  to  the  above  procedure  are 


t  1    (k)   (k)  m   L  ^3  (k)   1  (k)^   .  "I 

Lc(-r   >  r  >  s)  =  LDi3\^2  5     +  2  £  y     J 

x  [D23(^s(k)-21£(k))-1]  >  C^'3-5) 
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V£(k)> 
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"Fi-.^J-wfil^r^ria 
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*•  *•"     „<i  /        /  ~i  +  ~i-l 


i    (k)x 

/13I  \S      '  =  exp 

12  3 J 


{-*  &['(W-.«*) 


fe(^ 


4v    ^ 


-M       n       £l2 


-J        J     *• 


+  Q  + 


-12 


>, 


2*  V    » 

NTT 


x    /  Si  +  Si-i 


-12 


]}■ 


(4.1.29) 


V(r)   =  kK 


[(?) 


12 


I,    K  =  14.04  x    10"        erg,    a  =  2.56  x   10"     cm  , 


(2.1.2) 


and 


T    (Q)d2S  =  — i—  exp 

^2*7 


x       \        1 


r 


27  '  /  JzH- 


exp 


Tt/ 


2/ 


idQxdQy 


(4.3-9) 


The  procedure  for  computing  the  approximation  g  (r,p;M)  to 

g  (r,0)  is  the  same  as  that  for  g  (r,0;M)  except  that  steps  3  and  4 
e  ±-cL  c  Ld 


are  modified. 
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3.   Generate  M  independent  random  vectors  Q  =  (pA,0A,O)  , 

~t        <t    H     t 

t  =  1,2, .. .  ,M,  according  to  a  probability  T   (Q)d  Q. 
h.      Form  the  average 


(r^M)  =  |   Z   A(r[k),  q£k),  ^ )  r1^)   .  (B.3) 


t=l 


The  additional  equations  needed  for  the  computation  of  g  (r  ;M)  are 


,    (k)        (k)     Ax        1  Tt    /    (k)        (k)     n,       T    (    (k)        (k)       As 
(r        ,    1      S   Q)   =  o     L   (£        >    4       >   Q)   +  L0 (£       *    4       *    -Q) 


r    /    (k)  (k)     ns        _    ,    (k)  (k)        _s  " 

+  Le(rv    ',    -qv    ',    Q)   +  Le(rv      ,    -qv    ',    -Q)        , 


(U.3.210 


l.(sW,  S0J^,  ■  2  [,13  (S,  £  3(k)  ♦  I  ^(_k)  -  "uCS.  2(k))  ]  > 


D13(Q,q(k!))    =   exp 


*  £[Tte*«+ *) 


;+    lfV 


\fSi 


s^Tjt 


+  v 


fe^H-OK^H] 


(^.3.21) 


W;.) 


(1^.3.22) 


and 


T   (Q)d2Q  =  — %  exp 
e   ~        ~        6ly 


Ct) 


%  d|Q  "  ? -S  6<J  '  C0S  SQ  ^  1 


(^•3.19) 
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Table  1  gives  the  results  of  the  calculations  of 
\~3g  (r   ;20,000)  and  \"3g  (r^^O^OOO)  along  with  the  values  of 
W  (r   )  computed  in  [6],   The  quantity  \~^g  (r   ;20,000)  = 
\~3g  (r12;20,000)  +  \"3g  (r12;20,000)  is  shown  graphically  as  a 
function  of  r  p  in  Figures  17  through  25.   The  lengths  of  the  error 
bars  in  these  figures  are  twice  the  sum  of  the  standard  deviations 
in  \"3g  (r12;20,000)  and  \"3g  (^120,000). 

Table  2  presents  the  values  for  the  third  virial  coefficient 
calculated  from  the  formula 


c         3 


J  g1(r12)  W2(ri2)  d3r2  -  Ubg2    .        (2.3.20) 


The  integration  over  r,p  was  performed  by  Simpson's  rule  using 
g-,  (r12;20,000)  for  g,(r-Ip)  and  taking  the  values  of  Wp(r12)  and 


bp  calculated  in  [6]. 
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Table  1.1.   NUMERICAL  VALUES  OF  \"3g  (r12;20,000) , 

k~3gQ(r -^-,20,000),   AND  W2(r12)  FOR  T  =  5°K. 


T 

12 
a 

^gc(r12;M) 

Std.    Dev. 

^3  ge(r12>'M) 

Std.    Dev. 

W2(rL2) 

1.0 

0.11161 

0.05+33 

0.01^93 

O.O5076 

O.Ul620 

1.3 

+0.00768 

0.05091 

0.02307 

0.03620 

1.39898 

1.5 

-0.09163 

0.0^704 

+0.0260+ 

0.05+86 

1.5252  + 

1.7 

-0.155^2 

0.0^17+ 

-0.00652 

O.O38U2 

1.+1138 

1.8 

-0.16863 

O.O3699 

-0.02636 

0.03781 

1.33505 

1.9 

-0.10887  , 

0.03502 

,   -0.03+88 

0.05512 

1.26+86 

2.1 

-0.10791  '■ 

0.03155 

+0.02816 

0.0+026 

1.15^33 

2.3 

+0.06+28 

0.02298 

-0.01*913 

0.039^0 

I.O8719 

2.5 

0.15227 

O.OI892 

-0.02889 

0.0^606 

1.0^950 

2.8 

0.20509 

0.01092 

+0.0U91+ 

0.0++30 

1.02273 

Ill 


Table  1.2.   NUMERICA1  VALUES  OF  k'^g   (r   ;20,000), 

\"3ge(r19;20,000),  AND  W2(r12)  FOR  T  =  10°K, 


r 
J12 

a 

3gc(r!2^M) 

Std,    Dev. 

75  ge(ri2'M) 
K 

Std.    Dev. 

W?(r12) 

0.6 

0.30204 

0.05718 

-0,10363 

O.O5672 

0.00000 

0.7 

0.43065 

0.05918 

0 . 00001 

0.8 

0.22653 

0.05944 

-O.H826 

0,04663 

0,00845 

0.9 

0.20^32 

0.05659 

0.11992 

1.0 

0.21626 

0.05623 

+0.02633 

0,03522 

0.43502 

1.1 

+0.00022 

0.05535 

0,84777 

1.2 

-0,09561 

0.05339 

0,11547 

0,03143 

1,17350 

1.3 

-0,13564 

0.05041 

0,00588 

0,03157 

1,33885 

1.4 

-0.26533 

0.04638 

0.02020 

0.03433 

1,38106 

1.5 

-0.31497 

0.04414 

0,03717 

0.03018 

1=34965 

1.6 

-O.43696 

0.03915 

0.00436 

0.03249 

1,28194 

1-7 

-0.^3063 

0.03585 

+0.03937 

0,06796 

1,21369 

1.8 

-C, 48012 

0,03213 

-0.12736 

O.O3236 

1,15856 

1.9 

-0.45266 

O.02916 

1.11587 

2.0 

-O.41747 

O.02496 

+0.0106l 

0,03334 

1 . 08438 

2.1 

-0,25859 

0.02162 

1,06163 

2.2 

-O.I878O 

0.01845 

+0.04063 

0,03334 

1.04551 

2.3 

-0.09946 

0.01625 

-0,01458 

0,03264 

1,03408 

2.4 

+0.00139 

0.01320 

+0.02194 

0,03164 

1.02588 

2.5 

0.07724 

0,01100 

+0,05166 

0.03082 

.1 ,  01991 

2.6 

0.10874 

0.00953 

-0,03222 

0.03318 

1.01550 

2.7 

0.l4l66 

0.00755 

+0.00224 

0,03226 

1.01220 

2.8 

0.13137 

0.00595 

+0.00217 

0.032.17 

1 . 00970 

2.9 

O.14632 

0,00471 

1,00778 

3.0 

0,11885 

O.OO367 

-0.00508 

0.033^8 

I.OO629 

3-2 

O.08919 

0.00214 

1.00420 

3.* 

0,0667k 

0.001^9 

1,00289 

3.6 

0.04890 

c, 00099 

1.00203 

3-8 

0.03304 

0,00072 

' 

1,00145 

4.0 

0,02155 

0.000^7 

1  .  '  0106 
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Table  1.3.   NUMERICAL  VALUES  OF  X   3g  (r   ;20,000), 

\"3ge(r12 ;20,000),  AND  W  (r^)  FOR  T  =  20°K. 


r!2 

0 

^Jgc(rl2;M) 

Std.  Dev. 

^3  ge(r12;M) 

Std.  Dev. 

W^r^) 

0.6 

1 . 7702 

0.08450 

-O.62796 

0.07301 

0.00000 

o.t 

1.5075 

0.08448 

-0.21953 

0.05l6l 

0.00000 

0.8 

1.4115 

0.08353 

-O.HO3O 

0.04157 

0.00827 

0.9 

1.1673 

0 . 08096 

-0.08044 

0.03704 

0.14649 

1.0 

O.78071 

O.07988 

+O.O2763 

0.03097 

0.52629 

l.i 

0.53W 

O.07558 

0.02948 

0.03080 

0.94670 

1.2 

+O.21673 

O.07107 

0.08262 

0.02909 

1.19770 

1.3 

-0.013^6 

0.06457 

O.O3672 

0.03007 

1.26763 

1.4 

-0. 31^30 

0.05738 

0.01662 

0.03105 

1.24238 

1.5 

-0.48981 

0.05180 

+0.03585 

0.02992 

1.18817 

1.6 

-0.76050 

0.04506 

-0.02124 

O.03089 

1.13647 

1.7 

-0.85461 

O.03851 

+0.05727 

0.02978 

1 . 09690 

1.8 

-0.90454 

O.03174 

-0.02649 

0.03009 

1.06884 

1.9 

-0.84328 

O.02588 

-O.O2634 

0.03156 

1.04937 

2.0 

-0.722.40 

0.02248 

-0.03736 

0.03103 

1.03591 

2.1 

-0.58432 

0.01759 

-O.06747 

0.03031 

1.02652 

2.2 

-0.38289 

0.01473 

-0.01376 

0.03141 

1.01986 

2.3 

-0,23783 

0.01188 

-0.01845 

O.O2985 

1.01508 

2.1+ 

-0.08979 

0.00993 

-0.01048 

0.03039 

1.01159 

2.5 

-0.01504 

0.00764 

+0.01846 

0.03114 

1.00901 

2.6 

+0.03631 

0.00593 

-0.04298 

0.03033 

I.OO708 

2.7 

0.06169 

0.00461 

-0.01156 

0.02973 

1.00561 

2.8 

0.06026 

O.OO369 

+0.03174 

0.02972 

1.00449 

2.9 

0.05736 

O.OO298 

+0.02549 

0.03214 

1.00362 

3.0 

0.05746 

0.00232 

-0.00994 

0.03077 

1.00294 

3-2 

0.04262 

0.00149 

I.00199 

3-h 

0.03200 

0.00101 

1.00137 

3.6 

0.02172 

0.00069 

1 . 00097 

3.8 

0.01541 

0.00051 

1.00070 

U-.O 

0.01030 

0.00038 

1.00051 
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Table   1.4.      NUMERICAL  VALUES  OF  A~3g    (r      ;20;000), 

V"3ge(r12;20,00°- ^   AND  W2^rl2^   P0E  T  =   3°°K' 


12 


\ 


^  gc(.ri2;M)        Stdc    Dev,    j    A  gJ 


v 


:.r1?;M) 


Std,    Dev, 


wp'ri?.) 


-4- 


0.6 

0.7 
0,8 

0.9 
1.0 

1.1 

1.2 

1  =  3 
1.4 

1.5 
1.6 

1.7 

1.8 

1.9 
2.0 
2.1 
2.2 
2.3 

2.5 
2,6 

2.7 
2.8 

2.9 

3.0 
3.2 
3.4 
3.6 
3.8 
4.0 


3 . 7089 
3.1127 

2.9934 

2.6:- 
2,0920 
1.6787 
1.0346 
+0.50982 
-O.I29O8 
-0.62747 
-1,0382 
-1.1272 
-1.2507 
-I.I69I 
-0,96546 
-0,71601 
-0.46636 
-0,28688 

-0,15784 

-0,07918 

-c, 03204 

+0,01384 

0,01602 

0,02666 

O.C1990 

0.01861 

0.01 

O.OC919 
0,00681 

0.00479 


0,12325 
O.II960 

0, 11860 

0.11824 

0.11073 
0,10490 
0,09486 
O.O869O 
0.07376 
0.06234 
0.05072 
0.04l8l 
0,03376 
0.02662 
0,02114 
0,01693 
0,01305 
0,01032 
0,00818 
0,00650 
0.00505 
O.OO396 
0.00310 
0.00248 
0.00198 
0.00131 
O.OOO89 
0.0006^ 

0.C0048 

Oo 00037 


-0.55485 
-0.33827 

-0,06397 
-O.08559 

-O.OL259 

+o.oo4o3 

-0.01300 

+o.o6i48 

0.01243 

0,01058 

+0,03607 

-0,01810 

-0,03943 

-0.00269 

+0,06968 

0.03424 

+0.00379 

-0.00871 

-0.00372 
-0,00138 

+0,0284.1 

-0,03363 
+0,034146 
-0,01367 
-0,03350 


0.07420 
0.05606 

0.04174 

0.03514 

0.03C71 

0.03031 
0.03004 
0.03104 
0.03083 
0,03122 
0.03135 
0,0295? 
0.0306l 
Oo03l43 

0,03079 
0.03105 
O.O3.II6 
0.03214 
0.03316 
0.03094 
O.O3255 
O.C3069 
0.03112 
0.0312.9 
0.0309  3 


0.00000 
0,00001 

0,01251 
0.17920 
0.58907 
0,98567 
1.18254 

1,21059 

1.17199 
1.12496 
1.08763 
l,06l4l 
1,04.354 

1.03135 

1 . 02292 

1  -  01701 
1.01280 
1 . OO976 
1,0075? 
1.00586 
1 . 00462 
I.OO367 
I.OO294 
1,002  38 
1 , 00193 
1.00131 
1.00090 
1.00064 

1.000^6 
1.00034 
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Table  1.5.   NUMERICAL  VALUES  OF  \"3g  (r12;20,000) , 

\"3ge(r12;20,000),  AND  Wg(r^)  FOR  T  =  40°K. 


12 
o 

Jj6c(r12;M) 

Std.    Dev. 

3ge(r12;M) 
K 

Std.    Dev. 

W2(r12) 

0.6 

6.5168 

O.I6779 

-I.34129 

0.09346 

0.00000 

o.T 

5.8136 

0.16632 

-0.30960 

0.05813 

0.00002 

0.8 

5.2420 

0.16122 

-0.12666 

0.04290 

0.01505 

0.9 

4.6279 

0.15673 

-O.O6758 

0.03598 

0.21818 

1.0 

3.7103 

0.14460 

+0.05650 

0.03208 

0.66868 

1.1 

2.8058 

0.13353 

0.05465 

0.03152 

1.03258 

1.2 

2.0697 

0.12226 

0.03208 

o. 03006 

1.16868 

1.3 

1 . 1047 

0.10556 

+0.02770 

0.03180 

1.16784 

1.4 

+0.10555 

O.O8817 

-O.O2226 

0.03205 

1.12805 

1.5 

-O.51087 

0.07306 

-0.00293 

0.03048 

1.09045 

1.6 

-1.0742 

0.05946 

-0.05644 

O.03306 

1.06301 

1.7 

-1.3525 

0.04512 

+0.03024 

O.03016 

1.04421 

1.8 

-1.4343 

0.03477 

0.01522 

0.03201 

1.03146 

1.9 

-1.3756 

0.02575 

0.01866 

O.03218 

1.02273 

2.0 

-1.1223 

0.02022 

O.OO65I 

0.03H3 

1.01667 

2.1 

-0.81214 

0.01576 

0.03297 

0,03173 

1.01241 

2.2 

-0.56541 

0.01234 

0.02112 

O.03182 

1.00936 

2,3 

-0.3^673 

0 . 00976 

+0.00832 

O.03198 

1.00715 

2.4 

-O.20666 

0.00754 

-0.03604 

O.03083 

1.00552 

2.5 

-0.12062 

0.00592 

+o.ooi46 

O.03188 

1,00431 

2.6 

-O.06254 

0.00459 

-0,04644 

O.03171 

1.00340 

2.7 

-O.03496 

0.00363 

+0.04956 

O.03278 

1.00270 

2.8 

-O.OH60 

0.00289 

0.00766 

0.03141 

1.00217 

2.9 

-0.00110 

0.00231 

+0.02028 

0.03253 

1.00176 

3.0 

-0.00212 

0.00189 

-0.06352 

O.03143 

1.00143 

3-2 

•1-0.00275 

0.00127 

1.00097 

3^ 

0.00410 

0.00087 

1.00067 

3.6 

0.00224 

0.00065 

1 . 00048 

3.8 

0.00191 

0.00048 

1.00034 

4.0 

0.00125 

O.OOO38 

1.00025 
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Table   1.6.      NUMERICAL  VALUES  OF  \"3g   (r10 ;20,000), 


c^12j 


X"3ge(r12;20'°00)>  AND  W2(rl2)   F0R  T  =   5°°K' 


12  1  /  M\ 

^  gc(r12;M) 


Std.    Dev. 


\ 


3  ge(r12'M) 


Std.    Dev. 


w 


0.6 
0.7 
0.8 

0.9 
l.o 

l.i 

1.2 

1.3 
1.4 

1.5 
1.6 

1.7 
1.8 

1.9 
2.0 
2.1 
2.2 

2.3 
2.1+ 

2.5 
2.6 

2.7 
2.8 

2.9 
3.0 
3-2 

3.4 

3.6 

3.8 
4.0 


9.566^ 

8.7736 
7.3618 
6.6846 

5.5319 
4.4098 
2.8499 
I.5876 
+O.33265 
-O.48695 
-1.1023 
-1.6007 
-1.6056 
-1.5481 
-I.2836 
-O.91036 
-0.60499 
-0.40800 
-0.24114 
-0.l4l69 
-0.09400 
-0.05507 
-O.O3636 
-0.02176 
-0.01912 
-0.00742 
-0.00353 
-O.OOI89 
-0. 00182 
-0.00144 


0.21727 
0.21338 
O.20678 
0.19733 
0.18309 
0.16844 
O.14929 
0.12432 
0.10413 
0.08494 
O.06671 
0.04780 

0.03575 
0.02594 
0.01960 
0.01526 
0.01184 
0.00940 
0.- 00717 
O.OO561 
0.00442 
0.00350 
0.00276 
0.00221 
0.00182 
0,00124 
O.OOO88 
O.OOO65 
0.00050 
0.00039 


-O.8389O 
-0.41851 
-0.22029 
-0.02030 
+0.02428 
0.01908 
0.02397 
0.00415 
+0,04323 
-O.OO670 
+0.01203 
-0.03959 
+0.03151 
0.00077 
+0.02970 
-O.OH85 
-0.00923 
-0.05462 
+0,02053 
+0.00280 
-O.OO663 
-0.01^96 
+0,03584 
+0,03949 
-0, 06284 


0.08907 
O.06036 
0.04423 
0.03489 
0.03138 
O.03063 
O.03038 
0.03142 
O.03178 
0,03169 
O.03228 
0,03172 
O.03119 
0.03222 
O.03230 
O.03153 
0,03129 
0,03129 
O.03092 
O.03235 

0.03239 
0,03231 

O.03251 
O.03107 
0.03284 


!  1 


0.00000 
0.00001 

0.01763 
0.24798 
0.71891 

1.05350 
1.15232 
1.13778 
1.10136 
1 . 07088 
1.04933 
1.03467 
1.02472 
1 , 01789 
1.01315 
00980 
1,00740 
1,00566 
1,00438 
1,00342 
1,00270 
1,00215 
1.00172 

1.00140 

i,oon4 
1,00077 

1.00053 

1 , 00038 

1.0002.7 
1.00020 
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Table   1.7.      NUMERICAL  VALUES  OF  \"3g   (r      ;20,000), 

X'^S^r^; 20, OOO),   ANDW2(r12)   FOR  T  =  75°K. 


r12 
a 

^gc(r12^M) 

Std.  Dev. 

^3  6e(r]2;M) 

Std.  Dev. 

Vr12' 

0.6 

0.00000 

0.7 

17.099 

0.35252 

-0.57189 

0.06499 

0.00005 

0.8 

14.592 

0.33201 

•- 0.08237 

0.04288 

0.03028 

0.9 

13.849 

0.32565 

-O.O3IO2 

0.03544 

0.34036 

1.0 

9.9851 

0.28127 

+0.01446 

0.03462 

0.82732 

1.1 

8.1720 

O.25669 

+0.07193 

0.03240 

1.07206 

1.2 

5 . 6260 

0.22191 

-0.02478 

0.03356 

1.11255 

1.3 

3.4497 

0.18889 

-O.05651 

0.033++ 

1.09178 

1.4 

+1.3873 

0.14718 

+0.02364 

0.03500 

1.06573 

1.5 

-O.29434 

O.IIO83 

0.05378 

0.03440 

1.04574 

1.6 

-I.2900 

0.08048 

+0.04425 

0.03376 

1.03187 

1.7 

-I.9094 

0.05446 

-0.01024 

0.03336 

1.Q2246 

1.8 

-2.0339 

0.03728 

-0.01130 

0.03456 

I.OI606 

1.9 

-1.8976 

0.02542 

+0.04435 

0.03305 

1.01166 

2.0 

-1.5165 

O.OI976 

-0.00107 

0,03400 

I.OO858 

2.1 

-1.0853 

0.01530 

+0.04414 

0.03454 

1 . 00641 

2.2 

-0,7H91 

O.OII58 

-0. 00974 

0.03395 

1.00485 

2.3 

-0.48689 

O.OO897 

+0.00430 

0.03431 

1.00371 

2.4 

-O.33031 

0.00688 

0.00171 

0.03486 

1.00288 

2.5 

-O.22195 

0.00540 

0.02422 

0.03323 

1.00225 

2.6 

-0.15578 

0.00422 

+O.OI697 

0.03309 

1.00178 

2.7 

-0.11421 

0.00340 

-0.00746 

0.03451 

1 . 00142 

2.8 

-0.07068 

0.00268 

-0.00444 

0.03284 

1.00114 

2.9 

-0.05547 

0.00218 

-0.04366 

0,03403 

1 . 00092 

3-0 

-0.04714 

0.00182 

-0,05734 

0.03413 

1.00075 

3-2 

-O.02699 

0.00124 

1.00051 

3.4 

-0.01723 

O.OOO89 

1.00035 

3.6 

-0.01146 

O.OOO69 

1.00025 

3.8 

-0.00782 

0.00051 

1.00018 

4.0 

-0.00605 

0.00041 

1.00013 
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Table   1.8.      NUMERICAL  VALUES  OF  X_:5g   (r;L2;20,000), 

\"3g   (r^;20,000),   AND  W^r^)   FOR  T  =   100°K. 


12 


Tgc(r12'M) 


Std.    Dev. 


K 


h  ge(r12'M) 


Std.    Dev. 


W^Pjg) 


0.6 

0.7 
0.8 

0.9 
1.0 

l.i 

1.2 

1.3 
1.4 

1.5 

1.6 

1.7 

1.8 

1.9 
2.0 
2.1 
2.2 

2.3 

2.4 

2.5 
2.6 

2.7 
2.8 

2.9 
3.0 
3.2 
3-4 
3.6 
3.8 
4.0 


28.866 
26.301 

23.154 

20.788 

16.211 

12.226 
8.7524 
4.9980 
2.5440 

+0.15175 

-1.2950 

-2.2048 

-2 . 4492 

-2.1632 

-1.6750 

-1.1995 

-0.80334 

-0.57866 

-0.40427 

-0.27668 

-0.19290 

-0.14304 

-O.IO312 

-0.08112 

-0.06210 

-0.04i63 

-0,02636 

-O.O.1764 

-0,012.12 

-0.00962 


0.52539 
0.50866 
0.47882 
0.45434 
0.39965 
0.35103 
0.30203 
O.2369O 
0.19197 
0.14091 
O.O9839 
0.05886 
O.O3686 
0.02597 
0.02022 
0.01549 
0.01158 
O.OO898 
O.OO698 

0.00539 
0.00426 
0.00341 
0.00272 

0.00223 
0.00182 
0.00128 
0.00092 
O.OOO69 
0.00054 
0.00043 


-1.42417 

-0.78737 
-0.18906 

-0.01739 
+0.04359 
+0.00389 
-0.02168 
+0.05336 
+0.01317 
-O.03295 
-0.02260 
-O.07566 
+0.0621.4 
-0.01889 
-0.01509 
-0.06320 
-0.06583 
+0.02067 
-0.00948 

+0.01531 
-0.00674 
-0.02822 
-O.OO883 
-0.03814 
+0.03840 


0.11788 
O.06773 
0.04i4l 
0.03523 
O.03273 
0.03494 
0.03340 
0.03437 
0.03412 

O.03369 
O.03360 
O.03561 
0.03393 
0.03417 
0.03374 
O.03469 
0.03448 
0.03428 
O.03366 

O.03375 
0.03428 
O.03388 
0.03550 
O.03436 
0  03423 


0.00000 

0.00005 

0.04034 

0.40579 
0.87784 
1.06646 
I.08760 
1 . 06906 

1.04907 
1.03^11 

1.02378 
1.01678 
1.0.12' 01 

1 . 00872 
1 , 00643 

1,00480 
1 . OO364 
1.00278 
1.00216 
1 , OOI69 

1.00133 
1.00106 

I.OOO85 
1 . OOO69 
1,00056 
1,00038 
1.00026 

.1 ,  00019 
1.00014 

1,00010 
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Table   1.9.      NUMERICAL  VALUES  OF  \"3g   (r      ;20,000), 

X"3ge(r12;20'00°)>  MD  W2^r12^   F0R  T  =  2T3-l8°K. 


12 


^gc(r12^M) 


Std.    Dev. 


36e(r12;M) 


Std.    Dev. 


W2(r12) 


0.6 

0.7 
0.8 

0.9 
1.0 
1.1 

1.2 

1.3 
1.4 

1.5 

1.6 

1.7 
1.8 

1.9 


2.8 

2,9 
3o0 
3.2 
3-4 
3.6 
3-8 
4.0 


125.44 

115.17 

96.324 
81.748 
61.584 
45.317 
31.259 
19.126 
9.6566 
+2 . 3783 
-2 . 0686 
-3.8415 
-3.7872 
-3.1047 
-2.2944 
-1.6418 
-1.1846 
-0.84689 
-O.60708 
-0,45287 
-0,35206 
-O.26398 
-0.20450 
-O.15736 
-0,12176 
-O.08107 
-0.05550 
-0. 03776 
-0.02727 
-0, 02110 


2.0011 

1.9319 
1.7268 

1.5693 

1.3248 

I.0947 

0.89040 

O.66217 

O.47061 

O.28913 

0 . 14422 

0.07019 

0.04234 

0.03182 

0.02359 

0.01739 

Oo01324 

0.01005 
0.00770 
O.OO605 
0.00483 
O.OO389 

0.00315 
0.00256 
0.00213 

0.00153 
0.00113 
O.OOO85 
0.00068 
O.OOO56 


-I.81087 
-O.58497 

-0.04505 
-0.04151 
-0.03488 
-O.04590 
+0.01854 
+O.06585 
-0.02573 
+0.07952 
+0.02492 
-0.04085 
+0.00200 
-0.02403 

+0.00917 
+0.03585 
-0.00987 
+0.11174 
-0.00914 
+0.01644 
O.05927 
0.00352 
+O.00066 
-0.02901 
+0.01181 


0.11496 
O.06149 
0.0i+28l 
O.03984 
O.04038 
0.04040 
0.04181 
0.04014 
0.04132 
0.04200 

0.03909 

0 . 04049 
0.04228 
0.04070 
0.04094 
0.04239 
0.04263 
0.04l4l 
0.03990 
0.04090 
0.04145 
0.04258 
0.04304 
0.04089 
0.04266 


0.00000 
0.00119 
O.18533 
O.72582 
O.98181 
1.03326 

1.03307 
1.02481 
1.01745 
1.0L213 
1.00848 
1.00600 
1.00430 
I.00313 
I.00231 
1.00173 
1.00131 
1.00101 
1.00078 
1.00061 
1.00048 
1 . OOO38 
1.00031 
1.00025 
1.00020 

1 . oooi4 

1.00010 

1.00007 
1.00005 

1.00004 
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Table  2.   NUMERICAL  VALUES  OF  THE  THIRD  VIRIAL 
COEFFICIENT  C(T)  AND  THE  CLUSTER 
INTEGRAL  b„ . 


T(°K) 

c  (  cm  ^ 

b2(cm3) 

\mole  J 

10 

515 

1,-68 

20 

286 

-k55 

30 

226 

-1.50 

ko 

195 

-^.15 

50 

178 

-7^27 

75 

162 

-16,62 

100 

150 

-27*58 

273.18 

107 

-13O087 
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FIGURE    17.        X'^  AND  SAMPLING  ERROR  FOR  T=5°K. 
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FIGURE    18.       X'3g     AND  SAMPLING  ERROR   FOR  T=10°K. 
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FIGURE    19.       X"3g1  AND  SAMPLING  ERROR   FOR  T=20°K. 
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FIGURE    20.      X"°g1  AND   SAMPLING   ERROR   FOR  T=30°K. 
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FIGURE    21.      X"3g1  AND  SAMPLING   ERROR  FOR  T=40°K. 
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FIGURE    22.       X"3g1  AND  SAMPLING  ERROR  FOR   T=50°K. 
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FIGURE   23.      \-3g!    AND    SAMPLING    ERROR   FOR   T  =  75°K 
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FIGURE    24.      X"^   AND   SAMPLING    ERROR    FOR    T=100°K 
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FIGURE  25.      \~*q±  AND   SAMPLING   ERROR   FOR   T=  273.18  °K 
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